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ABSTRACT 

This  thesis  presents  two  methods  of  ship  classification 
with  IR  images.  The  first  method  is  the  Fourier  Coefficient 
method  which  transform  the  sample  points  of  a  superstructure 
profile  of  a  ship  to  the  spatial  frequency  components.  The 
shape  of  the  coefficient  curve  can  be  used  to  classify  the 
type  of  a  ship  from  8  different  categories.  But,  the 
differences  is  so  minor  that  it  is  difficult  to  implement  a 
computer  program  to  recognize  it.  The  second  method  is  a 
B-spline  Coefficient  method  which  uses  the  uneven  spaced 
spline  coefficients  to  find  the  beginning,  the  peak,  and  the 
area  of  the  lumps  of  a  ship  for  classification.  This  method 
is  better  than  the  Fourier  Coefficient  method.  The  study  of 
These  methods  is  presented  here. 
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I .  INTRODUCTION 

It  is  quite  important  to  classify  the  ships  according  to 
their  types.  Classification  can  be  done  in  a  number  of 
different  ways.  The  simplest  is  the  visual  method  that  is 
prone  to  error.  Other  methods  for  classification  are  the 
Fourier  Coefficient  method  and  the  B-spline  Coefficient 
method.  In  these  methods  the  necessary  information  for  clas- 
sification can  be  obtained  from  the  superstructure  profile. 

The  Fourier  Coefficient  method  samples  the  superstruc- 
ture profile  at  every  chosen  points.  The  function  values  at 
the  sampling  points  are  transformed  into  the  spatial  compo- 
nents. The  logarithm  of  the  magnitude  of  these  components 
is  plotted  and  compared  with  the  standard  plot  to  recognize 
the  type  of  the  ship. 

In  the  B-spline  Coefficient  method,  the  spline  coeffi- 
cients along  the  X  axis  and  the  Y  axis  are  used  to  recon- 
struct the  superstructure  profile.  The  shape  of  the  curve  of 
the  spline  coefficients  is  in  some  way  similar  to  the  shape 
(the  position  of  the  lumps)  of  the  superstructure  profile. 
The  ship  classification  may  be  achieved  by  recognizing  the 
beginning,  the  peak,  and  the  area  of  the  lumps. 
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II .  PREPROCESSING 

Preprocessing  is  the  procedure  to  obtain  the  superstruc- 
ture profile  of  a  ship  from  the  IR  image.  Then,  Fourier 
coefficient  or  spline  coefficient  methods  can  be  used.  The 
details  of  the  preprocessing  procedures  are  a  follows. 

A.   DATA  COLLECTION 

The  data  consists  of  the  IR  image  of  eight  different 
types  of  ships. 

1.  DD  -  Destroyer;  "HALL"  class. 

2.  Container;  The  class  is  unidentified. 

3.  Freighter;  The  class  is  unidentified. 

4.  AOR  -  Replenishment  oiler;  "WICHITA"  class. 

5.  LST  -  Tank  landing  ships;  "NEWPORT"  class. 

6.  FF  -  Frigate;  "GARCIA"  class. 

7.  CGN  -  Guided   missile   cruiser   (Nuclear  propulsion); 
"BAINBRIDGE"  class. 

8.  DDG  -  Guided   missile   destroyer;   "CHARLES  F.  ADAMS" 
class . 

These  images  are  taken  from  an  aircraft  which  is  flown 
at  a  height  of  500  feet  with  a  speed  of  400  knots  toward  the 
side  of  the  ship.  The  aspect  angles  for  these  images  are  90 
degrees  which  may  be  slightly  off  in  some  images.  The  inac- 
curacy of  the  aspect  angle  arises  from  the  fact  that  the 
photos  are  taken  while  the  aeroplane  keeps  on  moving.  All 
data  of  the  images  are  stored  in  a  digital  magnetic  tape 
with  256  by  64  bytes  per  image.  Thus,  the  number  of  bytes 
required  for  each  image  is  16384  and  each  byte  represents 
the  intensity  of  a  pixel.  For  each  record  of  the  image,  a 
label  is  coded  in  the  last  8  bytes  as  follows: 

1.   Byte  (16377)  =  Run  number  in  each  flight  which  passes 
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the  ship. 

2.  Byte  (16378)  =  Video  tape  time  code  when  the  data  is 
taken;  in  minutes. 

3.  Byte   (16381)  =   Video   tape   time   code;  in  seconds. 

4.  Byte  (16382)  =  Video  tape  time  code;  in  thirtieth  of 
a  second. 

5.  Byte  (16379)  =  Range  in  kilo-feet  which  is  the 
distance  measured  from  the  radar  it  may  have  an 
error  1  to  2  kilo-feet. 

6.  Byte  (16380)  =  Aspect  angle;  degrees  from  the  bow 
of  the  ship. 

7.  Byte  (16383)  =  Ship  class. 

8.  Byte  (16384)  =  ID,  1  =  for  training,  2,  3,  4,  5  =  for 
testing. 

The  run  number  and  the  time  code  together  uniquely 
define  a  specific  image  that  represents  a  single  TV  frame 
with  no  averaging.  In  addition,  the  time  interval  between 
the  end  of  one  image  to  the  beginning  of  the  other  is 
approximately  1.5  seconds.  Also,  there  are  inherent  random 
noise  in  the  record  which  arises  from  the  photo  instrument 
and  the  process  of  storing  them  on  to  the  digital  magnetic 
tape . 

B.   SOBEL  OPERATOR 

The  Sobel  Operator  technique  is  used  to  find  the  edge. 
To  determine  the  edge,  the  Sobel  Operator  uses  the  differ- 
ence of  gray  levels  of  the  pixels  in  a  3  by  3  matrix  as 
shown  in  Figure  2.1. 

a, b, c, d, e, f , g, h,  and  i  are  the  values  of  the  gray 
levels  at  the  position  of  (x-l,y),  (x,y-l),  (x-l,y),  (x,y), 
(x+l,y),  (x-l,y+l),  (x,y+l),  and  (x+l,y+l)  respectively.  The 
Laplacian  estimate  is  defined  as 

51  c^  f(x,yM(x+|,y)-[f(x+i,y)-f(x-f2,y)] 

5x2  {2-1} 
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The   basis   vector   for  all   directions   are   (a-2b+c), 
(g-2h+i),  (a-2d+g),  and  (c-2f+i).   Furthermore,   each  of  the 
basis  vector  is  convolved  with  the  image  as  follows: 
along  the  x  axis 


d*  =  [f (x-  l,y- 1)  +  2f(x,y-  |)+f(x+  |,y-Q] . 


-[f(x-l,y4-l)  +  2f(x;y+i)+f(x+!;y+i) 


(2.2) 


along  the  y  axis 


ty  =  ff(x+i,y-i)+2f(x+U)+f(x+i,y+i)" 


-[f(xH,y-!)+2f(xH,y)+f(x-!,y+l) 


(2.3) 


Since  the  magnitude  of  the  resulting  vector  is  the  abso- 
lute value  of  the  convolved  results,  the  edge  magnitude 
S(x,y)  [Ref.  1], 


S(x,y)  =  (d/+d,T 


(2.4) 


Note  that  the  Sobel  operator  does  not  use  the  gray  level 
at  the  position  (x,y).  The  advantage  of  using  a  Sobel  oper- 
ator over  others  is  that  the  resulting  edge  is  smoother  due 
to  a  3  by  3  matrix  approach.  If  we  compare  the  Sobel  oper- 
ator with  the  Laplacian  operator,  it  is  seen  that  the  Sobel 
operator  using  the  four  basis  vector  as  shown  above  will 
provide  more  accurate  reading  because  of  noise  reduction  in 
the  original  image.  Hence,  The  Sobel  operator  is  often  used 
in  the  preprocessing  operation. 
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b 

c 
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f 

g 

h 

i 

Figure  2.1    Sobel  Operator. 

C.   THE  USE  OF  A  SOBEL  OPERATOR 

When  a  Sobel  operator  is  used  at  the  edge  of  the  image 
frame,  the  pixel  level  which  lies  out  of  the  frame  will  be 
set  equal  to  that  of  the  adjacent  pixel  within  the  frame. 
The  original  image  of  a  guided  missile  cruiser  is  shown  in 
Figure  2.2.  Since  the  result  of  the  Sobel  operator  is  numer- 
ically greater  than  8  bits  range,  we  have  to  rescale  the 
result  back  to  8  bits  range.  This  is  achieved  by  deter- 
mining the  maximum  and  the  minimum  of  the  gray  level .  They 
are  then  used  to  rescale  the  gray  level  in  the  results  of 
the  Sobel  operator  as  shown  in  Figure  2.3.  In  some 
instances,  the  original  image  is  very  poor  as  shown  in 
Figure  2.4.  Attempts  to  determine  the  edge  of  this  image 
failed  as  shown  in  Figure  2.5. 


Figure  2.2    A  CGN  at  a  Range  of  45000  feet 
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Figure  2.3    Image  from  a  Sobel  Operator 


Figure  2.4    Blured  Image  of  a  LST  at  a  Range  67000  feet 


.      .-■•    ,' — -,-.-    ?■■■■  ■-:-;:-;'-  ■>-;:.?*:■■■.■:  -■    V!-..-  -v  ■'.-■..-'  ..  ■,„:'".  '■';.    -':..?'  .-;..VC    .-""     ■•.-■  ';>~  -••  .■     A>^!£3fVsfe?Srti3ii4SJ»l**->v-e'."-i; 

8^!?^'  "   '  Wfcau^j  i^4    ii^^^^^^^^ 

lB|||BlMG^3IIBfiiS3WBWHn^^                                                 aSBSSBtMKBi&ttiKiu»iiSiilSv^<^^^x9^Ki^^!SiS^!3Si 

Figure  2.5    Noisy  Image  from  a  Sobel  Operator 
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D. 


EDGE  THRESHOLD  STRATEGIES 


Edge  threshold  strategies   are  used  to  extract  the  edge 

profiles  from  the  Sobel  results.   In  this  case,   we  use  only 

the  pronounced   value  of  an   edge  element   at  x  if  g(x)   is 
greater  than  certain  threshold  value  [Ref.  2]. 


'    =<3(x,y)  if  g(xy);>,  threshold 


G(*,y)    { 


=  0 


otherwi  se 


(2.5) 


To  increase   the  contrast  of   the  image  to   a  silhouette 
form,  G(x,y)  is  defined  as 


f  =255 


G(*,y)    I 


=  0 


if    9(*,y)  ^threshold 


other  w  ije 


(2.6) 


The  choice  of  the  threshold  value  is  based  upon  the 
histogram  of  the  edge  image  as  shown  in  Figure  2.6.  The 
estimated  critical  gray  level  is  chosen  so  that  a  majority 
number  of  the  pixels  with  value  between  0  to  255  will  fall 
below  the  critical  value.  Alternatively,  histogram  equaliza- 
tion may  be  used  to  determine  the  desired  threshold  level  as 
shown  in  Figure  2.7.  In  this  case,  trial  and  error  method 
was  used  in  conjunction  with  the  above  method  to  obtain  the 
threshold  so  that  the  correct  profile  is  ascertained. 
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Figure  2.6    Histogram  of  Figure  2.3 
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Figure  2 . 7 


Cumulative  Distribution  of  the  Histogram 
in  Figure  2.6. 
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E.   HOW  TO  EXTRACT  THE  PROFILE 

The  edge  image  of  the  guided  missile  cruiser  shows 
little  variations  in  the  gray  level  as  shown  in  Figure  2.3. 
These  variations  are  caused  by  the  noise  in  the  original 
image.  In  this  case,  the  choice  of  the  threshold  value  is 
based  upon  both  the  histogram  and  the  cumulative  distribu- 
tion of  the  edge  image  so  that  it  contains  90  %  of  the 
pixels.  Therefor,  the  chosen  gray  level  is  110  and  the 
result  is  shown  in  Figure  2.8. 


Figure  2.8    Silhouette  of  a  CGN  in  Figure  2.3 


F.   HOW  TO  OBTAIN  THE  SUPERSTRUCTURE 

The  original  image  .  of  the  ship  is  taken  from  the  aero- 
plane with  different  displacement  from  waterline  to  the 
superstructure  in  a  rough  sea,  so  that  the  profile  of  the 
ship  with  respect  to  the  sea  surface  varies.  Thus,  we  have 
to  eliminate  some  information  in  the  image  of  the  ship  by 
considering  the  superstructure  only.  In  considering  the 
overall  ship  structure,  it  is  obvious  that  the  largest 
distance  is  between  the  bow  and  stern  span.  Therefore,  the 
bow  and  stern  points  are  located  first  in  the  program  as 
shown  in  Appendix  A.  Then  we  consider  the  slope  of  the  bow 
and  stern  of  the  ship  and  set   all  the  gray  values  below  the 
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slope  line  to  zero.  Hence,  it  results  in  a  superstructure  of 
the  ship  as  shown  in  Figure  2.9.  However,  in  some  images, 
there  is  a  lot  of  noise  in  the  background  which  cause  diffi- 
culty in  locating  the  bow  and  stern.  Under  these  circum- 
stances, the  dimensions  of  a  ship  are  estimated  by  trial  and 
error  method.  Then  the  gray  level  which  lies  outside  is  set 
to  zero  and  an  estimate  of  the  bow  and  stern  slope  can  be 
made . 


Figure  2.9    Superstructure  Profile  of  Figure  2.8 


G.   CONTOUR  FOLLOWING 

The  noise  in  the  background  of  the  original  image, 
yields  wide  edge  structure.  In  considering  this  factor,  the 
contour  following  procedure  tracking  the  inner  part  of  the 
image  in  Figure  2.9  gives  the  superstructure  profile.  The 
objective  of  this  contour  following  is  to  describe  the  bow 
and  stern  points  of  the  ship,  the  direction,  and  the  posi- 
tion of  the  edge  of  the  superstructure.  The  contour  tracing 
is  done  in  the  counter-clockwise  direction  which  compares 
pixel  value  of  0  or  255  in  a  3  by  3  matrix  in  the  following 
manner. 
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Starting  from  the  leftmost  point  in  the  superstructure 
image  in  Figure  2.9  the  contour  profile  tracing  is  accom- 
plished by  examining  the  neighbors  of  a  3  by  3  kernel 
located  at  the  curser  position  as  shown  in  Figure  2.10.  The 
curser  is  moved  along  the  profile.  All  successive  positions 
of  the  curser  constitute  the  contour  profile  of  the  ship. 
The  testing  procedure  is  explained  below. 

1.  Initialize  the  curser  position  to  the  beginning  of 
the  thresholded  image  with  the  gray  level  of  255  and 
the  estimate  direction  of  the  next  position. 

2.  Check  for  reaching  the  end  of  the  profile,  if  it  is 
at  the  ending  of  the  profile  then  stop,  if  not  go 
to  3  . 

3.  Check  for  the  estimate  to  see  whether  it  is  in  the 
direction  of,  North,  East,  South,  or  West.  If  the 
direction  it  is  North  then  go  to  4,  if  it  is  East 
then  go  to  5 ,  if  it  is  South  then  go  to  6,  and  if  it 
is  West  then  go  to  7. 

4.  Determine   the  next   position   with  the  gray  value  of 
255  in  the   counter-clockwise   direction   as  shown  in 
Figure  2.11.  Store  the  position  found  and  move  the 
curser   to   that   position.    Update    the    estimate 
direction  to  the  one  last  found;  then  go  to  2. 

5.  The  procedure  is  the  same  as  that  in  step  4  except 
search  pattern  is  shown  in  Figure  2.12. 

6.  The  procedure  is  the  same  as  that  in  step  4  except 
search  pattern  is  shown  in  Figure  2.13. 

7.  The  procedure  is  the  same  as  that  in  step  4  except 
search  pattern  is  shown  in  Figure  2.14. 

The  flow  chart  of  the  procedure  is  shown  in  Figure 
2.15,  and  the  detail  of  each  procedure  are  included  in 
Appendix  B.  The  testing  procedure  have  to  be  performed  in 
such  a  maner  that  the  resulting  contour  is  a  good  represen- 
tation of  the  superstructure  line.  The  result  of  the  contour 
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image  is  shown  in  Figure  2.16.  If  the  superstructure  in 
Figure  2.9  has  wide  edge,  we  can  not  find  the  contour  image. 
We  have  to  use  the  additional  step  which  is  called  the 
"Closing"  operation. 


Figure  2.10    The  3  by  3  Kernel 


Figure  2.11    The  Step  Checks  in  the  North  Direction 
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Figure  2.12    The  Step  Checks  in  the  East  Direction 


Figure  2.13    The  Step  Checks  in  the  South  Direction. 
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Figure  2.14    The  Step  Checks  in  the  West  Direction 
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Figure  2.16    Contour  Image  of  a  CGN  at 
a  Range  of  45000  feet. 


H.   CLOSING  OPERATION 

Some  results  from  the  superstructure  extraction  process 
are  discontinuous  because  the  gray  level  of  those  areas  of 
the  structure  is  less  than  the  threshold  value.  If  we 
decrease  the  threshold  value,  the  details  of  the  superstruc- 
ture are  effected.  It  is  necessary  to  use  the  "Closing" 
operation  which  consists  of  the  "dilation"  process  to  smooth 
the  superstructure  profile  used  the  "Erosion"  process.  All 
direction  are  dilated  similarly  which  causes  an  smoothing 
effect  on  the  edges,  The  superstructure  increases  in  total 
area.  Then,  use  the  "Erosion"  process  to  shrink  (subtract) 
the  dilated  part  in  all  directions,  thus  obtain  the  smoothed 
superstructure  with  the  appropriate  size.  The  Closing 
process  employs  the  dilation  process  which  yields  an  output 
image  from  an  input  image.  The  Dilation  results  is  shown  in 
Figure  2.17. 
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Figure  2.17    The  Process  of  Dilation. 

We  examine  the  gray  value  of  each  pixel  in  the  original 
image.  Begining  from  the  pixel  at  the  first  row  and  the 
first  column.   The  procedures  are  the  following. 

1.  Considering  one  pixel  in  the  original  image  with  the 
kernel  (B)  centered  there.  If  at  least  one  pixel  in 
the  kernel  has  a  value  of  255,  we  let  the  gray  level 
of  the  output  image  at  the  center  of  the  kernel  be 
255. 

2.  We  shift  the  center  of  the  kernel  1  column  to  the 
right.  Then  following  the  same  procedure  as  in  step  1 
until  the  last  column  is  reached. 

3.  We  shift  the  position  of  the  kernel  to  the  next  row 
and  starting  from  the  first  column.  Then,  following 
the  same  procedure  as  in  step  1  and  step  2  until  the 
last  row  and  the  last  column  is  reached. 

The  result  obtained  from  the  dilation  process  is  an 
image  with  enlarged  structure.  The  second  procedure  in  the 
Closing  operation  is  the  Erosion  process.  The  Erosion 
process  perform  the  same  procedures  as  the  dilation  process 
except  for  step  1.  If  every  pixel  in  the  kernel  are  255,  we 
let  the  gray  level  in  the  new  image  at  the  center  of  the 
kernel  be  255.  Otherwise,  it  will  be  0.  The  output  image 
obtained  will  have   smooth  edge  with  minor   change  occurring 
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in  the  edge  detail  as  shown  in  Figure  2.18.  In  this  case,  we 
use  an  kernel  of  size  3  by  3.  If  we  increase  the  size  of  the 
kernel,  the  details  of  the  image  are  decreased. 


Figure  2.18    The  Profile  after  Dilation  and  Erosion 

(Closing  Process). 


I.   PROFILE  ROTATION 

Often  in  the  contour  image,  the  first  and  the  last  point 
of  the  superstructure  are  not   at  the  same  horizontal  level. 
We   have  to   rotate  the   contour   image  by   setting  the   two 
points  to  the  same  horizontal  level.    How  to  rotate  it  from 
the  Y,  X  axis  to  the  Y',  X'  axis  is  shown  in  Figure  2.19. 

If  the  angle  value  Q   is  positive,  it  will  be 


X' 

Y' 


cosP      -sinfl 
8        cos£J 


sin 


If  the  angle  value  Q   is  negative,  it  will  b< 


X1' 

y 


cosQ        sinff 
n0       cos9. 


-si 


For  some  of  the  contour  profile,  part  of  the  profile 
after  rotation  will  be  out  of  the  image  frame.  Then  we  have 
to  shift  this  contour  down  by  20  pixels  position.  The 
rotated  profile  is  shown  in  Figure  2.20. 
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Figure  2.19    Rotation  Process 


Figure  2.20    Rotated  Profile  of  a  CGN  at 
a  Range  of  45000  feet. 


30 


Ill .  FOURIER  COEFFICIENT  METHOD 

We  have  obtained  the  ship  profile  in  the  previous 
chapter.  To  extract  features  out  of  this  profile  for  clas- 
sification purposes,  we  will  use  the  Fourier  Transform 
method. 

The  Fourier  transform  of  the  ship  profile  showed  that 
the  transform  coefficients  depend  upon  the  ship's  dimen- 
sions, its  superstructure,  and  the  distance  between  the 
camera  and  the  ship.  If  the  profile  f(x)  is  a  discrete 
function  with  128  sample  points.  The  discrete  Fourier 
transform  can  be  written  as 

M-/ 

F(u)  =l£f(x)cxp(-j27Tux/N)  (3>1) 

For  u,  x  =  0,  1,  2,  ,  N-l.  N  is  the  total  number  of 

samples . 

If  the  direct  calculation  of  the  discrete  Fourier  trans- 
form is  chosen,  the  number  of  complex  multiplication  and 
addition  will  be  equal  to  N  ;  i.e.  to  obtain  F(0)  would 
require  complex  multiplication  and  addition  N  times.  In 
order  to  reduce  the  computation,  we  use  the  fast  Fourier 
transform  algorithm.  Thus,  equation  3.1  [Ref.  2]  can  be 
separated  into  F even    (u)  and  F odd   (u) 

x=o 

F„,>)=iyf(2x+i)w;;x         <3.3> 
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W  =exp(-j27T/M) 


M  =  N. 
2 


(3.4) 
(3.5) 


F(U)=i[^e>)  +  ^>)^] 


(3.6) 


F(u+M)=l[Fr^»-F0£/,(u)^J 


(3.7) 


Using  this  method,  we  have  reduced  the  total  number  of 
complex  multiplication  and  addition  to  Nlog  N.  In  this  case, 
we  have  N  =  128,  thus,  the  total  number  of  complex  multipli- 
cation and  addition  will  be  896. 

First,  we  divide  the  rotated  profile  image  into  128 
divisions.  Since  the  distance  of  each  division  is  equal  to 
the  amount  of  the  pixel  between  the  bow  and  the  stern  of  the 
ship  divided  by  128,  we  use  the  distance  perpendicular  from 
the  horizontal  line  between  bow  and  stern  to  the  highest 
point  of  the  superstructure  as  the  sampled  values.  The 
result  of  the  fast  Fourier  transform  is  a  complex  number. 
Then  we  use 


M(u)  =  loq[|  +  G(u)] 


(3.8) 


G(u)  is  the  magnitude  of  F(u).  A  value  of  1  is  add  to  the 
magnitude  to  avoid  negative  logarithm  result,  the  results 
obtained  are  as  shown  in  Table  I  and 

1.  Figure  3.1  -  Destroyer 

2.  Figure  3.2  -  Container 

3.  Figure  3.3  -  Freighter 

4.  Figure  3.4  -  Replenishment  oiler 

5.  Figure  3.5  -  Tank  landing  ship 
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6.  Figure  3.6  -  Frigate 

7.  Figure  3.7  -  Guided  missile  cruiser 

8.  Figure  3.8  -  Guided  missile  destroyer 

On  the  results  minor  difference  of  the  shape  of  the 
Fourier  coefficients  can  be  noticed  visually.  But  it  is 
difficult  to  implement  a  program  to  detect  these  minor 
difference  in  shape.  Therefore,  a  Second  method  was  used  to 
handle  this  problem. 
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TABLE  I 
Logarithmic  Magnitude  of  Ships 


F(N) 

DD 

Cont. 

Frig. 

AOR 

LST 

FF 

CCN 

DDC   | 

(1E-1) 

(1E-1) 

(1E-1) 

(1E-1) 

(1E-1) 

(1E-1) 

(1E-1) 

(1E-1) 

(1E-1) | 

DC 

1.E1 

1.E1 

1.E1 

1.E1 

1.E1 

1.E1 

1.E1 

1.E1  | 

F(l) 

6.74 

7.23 

3.66 

4.28 

6.81 

6.97 

7.78 

8.  10  | 

F(2) 

2.  17 

7.18 

5.  17 

8.04 

2.64 

5.67 

6.37 

5.04  | 

F(3) 

3.21 

6.45 

7.88 

3.84 

3.35 

5.64 

2.38 

5.40  | 

F(4) 

3.17 

6.18 

3.76 

5.28 

1.28 

4.43 

7.6E-1 

3.26  | 

F(5) 

2.20 

5.09 

5.31 

2.45 

2.49 

2.  40 

3.53 

4.81  | 

F(6) 

1.88 

4.57 

3.  17 

2.66 

1.36 

1.52 

3.48 

6.02 

F(7) 

1.10 

2.03 

3.90 

2.28 

3.6E-1 

8. 1E-1 

3.44 

4.50  | 

F(3) 

1.50 

2.65 

4.58 

1.07 

1.74 

1.76 

3.35 

3.  14  | 

F(9) 

2.07 

1.83 

1.67 

1.32 

2.14 

9.9E-1 

3.40 

2.46  | 

F(10) 

1.  13 

1.57 

2.07 

3.39E1 

1.  19 

1.  13 

3.24 

1.05  | 

F(ll) 

4.  6E-1 

6.0E-1 

4.53 

1.03 

2.46 

1.  14 

2.62 

1.23 

F(12) 

1.6E-1 

5.6E-1 

2.23 

1.07 

1.73 

9.8E-1 

1.  80 

2.09  | 

F(13) 

5.4E-1 

1.62 

2.09 

5.0E-1 

1.78 

3.5E-1 

3.5E-1 

1.23 

F(14) 

9.  6E-1 

1.48 

2.09 

9.8E-1 

1.37 

7.4E-1 

6.7E-1 

1.62 

F(15) 

4. 1E-1 

1.25 

1.07 

7.CE-1 

1.14 

4. 1E-1 

1.19 

1.38  | 

F(16) 

2.9E-1 

1.77 

1.30 

1.05 

7. 1E-1 

1.36 

1.25 

1.30  | 

DD  =  Destroyer  at  range  77000  feet. 

Cont  =  Cointainer  at  range  28000  feet. 

Frig=  Freighter  at  range  40000  feet. 

AOR  =  Replenishment  oiler  at  range  78000  feet. 

LST  =  Tank  landing  ship  at  range  51000  feet. 

FF  =  Frigate  at  range  49000  feet. 

CCN  =  Guided  missile  cruiser  at  range  45000  feet. 

DDG  =  Guided  missile  drestroyer  at  range  41000  feet 
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A.   VARIATION  OF  SHIP  SUPERSTRUCTURE  WITH  RANGE 

One  of  the  practical  problem  of  using  the  ship  profile 
is  that  it  is  sensitive  to  range  variations.  Close  ship 
profile  has  more  details  than  the  far  away  ship  profile. 
The  dependency  of  the  geometric  size  of  the  profile  with  the 
range  is  discussed  in  this  section. 

Assume  that  the   object  is  centered  on   the  camera  axis. 
The  field  of  view  of  the  camera,   the  number  of  the  pixel  in 
the  image,  the  size  of  the  image,  and  the  size  of  the  object 
are  known.   Our  problem  is  to  determine  the  distance  between 
the  camera  and  the  object. 

1 .   Background 

The  camera  system  is  similar  to  a  human  eyes  system. 
The  light  reflected  from  the  object  goes  into  the  eyes.  The 
image  of  the  object  falls  on  the  retina,  the  signal  is  sent 
to  the  brain  in  some  electrical  from,  and  the  brain  changes 
it  to  a  from  that  human  can  perceive.  But,  in  the  camera 
system  film  or  senser  are  used  to  pick  up  the  image.  The 
function  of  a  camera  is  shown  in  Figure  3.9 


image 

\^^    lens     ^^^ 

object 

Figure  3.9    The  System  of  the  Camera 
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Figure  3.10    Range  Calculation. 

The  distance  between  the  lens  and  the  image  plane 
can  be  adjusted  in  order  to  have  a  clear  image  on  the  film. 
The  camera  has  a  field  of  view  (FOV)  angle  as  shown  in 
Figure  3.9.  The  object  has  to  be  in  the  field  of  view  of 
the  camera.  The  determination  of  the  distance  is  shown  in 
Figure  3.10 

For  simplicity  of  calculation  the  inside  of  the 
camera  is  flipped  to  the  same  side  as  the  object  as  shown  in 
Figure  3.10.  when  the  angle  of  the  FOV  is (2,  1/2  is  the 
half  of  the  full  image  size,  D/2  is  the  half  of  the  dimen- 
sion of  the  object,  X  is  the  distance  from  the  lens  to  the 
image,  and  R  is  the  distance  from  the  lens  to  the  object. 

Assume  that  the  distance  1/2,  distance  Im/2,  and 
angleQ^/2  are  known.  Then,  the  distance  R  can  be  determined 
by 


m_tana 
2.  Z 


(3.9) 


tana  =IL 
2      2X 


(3.10) 
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R  =_D  J_  (3.ii) 

2  tana 
2 

Assume  that  the  angle  resolution  of  the  pixel  of  the 
field  of  view  of  the  camera  is  equal  to  0.2E-3  radian  per 
pixel.  The  number  of  pixel  of  the  frame  is  equal  to  256. 
The  size  of  an  image  is  I  pixels.  The  dimension  of  the 
object  D  in  feet  is  known.  The  field  of  view  in  angle  is 


<2=(0.2E- 3)256 


(3.12) 


2 

where  d  is  in  unit  of  pixel. 


tana' =  I    d 


X  =  256  __d 

2     tanCL  (3.13) 


2      2   X  <3-14) 


tana=_I_  tancy 

2    256        2  (3.15) 


R=_Q_ 


2  tana;  (3.16) 

2 


R=128D 

ltanO.0256  (3.17) 


When  the  length  D  of  the  object  is  known,  from  equa- 
tion 3.17  we  can  estimate  the  distance  from  lens  to  the 
object  and  is  shown  in  Table  II 
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TABLE  II 
Range  Estimation 


|  CLASS 

I  ( pixel ) 

D  (ft) 

r'  (kft) 

R(kft) 

(R-r')IOO 

MEASURE 

KNOWN 

CALCULATE 

RADAR 

DIS 

R 

|  DD1 

96 

413 

21.766 

77 

71.71 

|  DD2 

80 

418 

26.  119 

85 

69.27 

|  A0R1 

107 

659 

30.787 

78 

60.53 

|  A0R2 

96 

659 

34.315 

85 

59.63 

|  LST1 

176 

522.3 

14.834 

51 

70.91 

|  LST2 

134 

522.3 

19.484 

57 

65.82 

|  CGN1 

147 

565 

19.213 

45 

57.31 

I  CGN2 

119 

565 

23.734 

55 

56.85 

|  DDG1 

126 

437 

17.337 

47 

63.  11 

|  DDG2 

90 

437 

24.247 

64 

62.08 

DD1.DD2  =  Destroyer  at  range  79000  and  83000  feet. 
A0R1,A0R2  =  Replenishment  oiler  at  range  78000  and  88000  feet. 
LST1,LST2  =  Tank  landing  ship  at  range  51000  and  62000  feet. 
CGN1.CCN2  =  Guided  missile  cruiser  at  range  45000  and  64000  feet. 
DDG1.DDG2  =  Guided  missile  destroyer  at  range  41000  and  64000  feet. 


The  error  distance  between  the  estimated  distance 
and  calculated  distance  in  percentage  is  ((R  -  R)/R)  100. 

2 .   Remark 

Calculated  distance  error  in  R  may  come  either  from 
the  pixel  measurement  in  an  image  or  from  the  angular  reso- 
lution estimation.  The  distance  that  is  stored  in  the  image 
label  has  an  error  of  about  1  to  2  kilo-feet.    If  the  angle 
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of  the  field  of  view  is  accurate,  we  can  estimate  the  number 
of  pixel  in  an  object  image.  Then  we  can  classify  the  object 
by  comparing  the  number  of  pixel  of  the  original  image  with 
that  of  the  test  image.  Some  known  system  parameters  can 
help  to  determine  the  range  of  the  target  ship.  The  problem 
is  that  existing  errors  in  the  system  parameter  causes  in 
the  larger  errors  in  the  estimated  range.  Consequently,  they 
are  not  very  useful  in  classifying  the  ships. 
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IV.  B-SPLINE  COEFFICIENT  METHOD 

Using  B-spline  coefficient  is  another  method  to  describe 
a  ship  profile.  This  method  uses  the  B-spline  coefficients 
to  determine  the  beginning,  peak,  and  area  of  lumps  which 
contain  significant  information  about  the  type  of  the  ships. 
The  comparisons  of  the  knot  position  (in  parametric  value) 
from  the  midships  to  the  peak  or  beginning  of  the  lump,  can 
be  very  helpful.  Different  ships  will  have  different  lump 
positions  and  sizes. 

A .   BACKGROUND 

A  spline  function  is  a  piecewise  polynomial  used  to 
interpolate  points.  This  kind  of  curve  is  smooth  and  the 
discontinuities  in  its  k-th  derivative  is  as  small  as 
possible.  In  this  case,  Cubic  spline  was  used,  where  the 
first  and  second  derivatives  for  any  set  of  interpolating 
points  are  continuous,  while  the  third  derivative  may  be 
discontinuous.  The  reason  for  using  Cubic  spline  is  to  keep 
the  third  derivative  discontinuity  as  small  as  possible  and 
the  curve  as  smooth  as  possible.  The  B-spline  calculation 
procedure  is  also  very  stable.  In  our  case  the  order  of 
spline  used  is  3,  while  the  1-st  and  2-nd  derivative  are 
continuous.  The  choice  of  4-th  order  spline  function  is  due 
to  the  fact  that  it  is  generally  sufficient  for  most  ship 
profile  curves.  Discontinuity,  in  a  sense,  may  be  stated  as 
the  jumps  of  the  third  derivative,  which  is  the  means  to 
control  smoothness  of  the  connecting  pieces. 


48 


B.   B-SPLINE  APPROXIMATION  WITH  FREE  KNOT 

As  mention  earlier,  the  B-spline  function  is  used  in  the 
spline  approximation  with  free  knot.  The  free  knots  is  some- 
times called  uneven  or  irregular  knots;  that  is  no  fixed 
number  of  knots  are  used.  Furthermore,  the  spline  position 
need  not  be  on  the  original  curve  so  as  to  minimize  the 
number  of  knot  position  while  preserving  most  of  the  detail's 
of  the  original  curve. 

First,  minimize  the  value  of  the  lack  of  the  smoothness 
77(c)  defined  as  [Ref.  3] 


n       n 


77(C)  =  Id  a-.C,.) 


j=l    i--k 


'•J 


(4.1) 


where  Ci   is  a  coefficient  of   spline  at  the   knot  position. 


a. ij     is  defined  as  follows 


-  uU) 


(*) 


a/;  =  M^+/(V+0)-M^+/(trO) 


(4.2) 


where  M 
as 


I ^+^    is  the  normalize   B-spline  function  and  defined 


M 


*+/, 


/,*+/(x)  =  (t/+4+rt/)Afr(t/,...It/+ik+/)qA(t;x) 


(4.3) 


and  G^(t;x),  t   are  defined  as  follows 


f=(t-x)*  =  (t-x)*   if  t^x 


'+' 


%M{ 


=  0 


if  t<x 


(4.4) 
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A^(Z;,Z,>/  ,...,Z/+£  )f(t)  stands  for  the  1-th  divided 
difference  of  the  function  f(t)  on  the  point  Z; , . . . , Z^ 
where  t  is  the  position  values  of  knots  in  term  of  Z  param- 
eter. Z  parameter  is  defined  as  follows 


Z(l)  =Z(I-I)+  [(X(l)-X(IH)f+(Y(l)-Y<IH))2]' 


(4.5) 


Z(0)  =  0 


(4.6) 


where   I    is   the    number   of    the   sampling    point   and 
1=1 , 2 ,  3 , . . . , m . 

Second,  the  smoothing  is  subjected  to  a  constraint 


6(c) 


<  5 


(4.7) 


where  S  is  the  smoothing  factor,  0(C)  i£  the  weighted  sum  of 
the  square  residuals  defined  as 

2 


5(c)=^.[y;~^c/m.a+/(x;)> 


(4.8) 


X: , Y-    are  the   values  of   X  and   Y   at  Z   parameter  of   the 
>ampl 
points  [Ref.  3]  defined  as 


sampling  points;   W-  is  a   weighting  factor  for  all  sampling 


W;  =  (5Yy) 


-2 


(4.9) 


where  W-  is  an  estimate  of  the  standard  deviation  of  Y- . 
Then,  the  value  of  S  is  in  the  range  m+V2m.  If  nothing  is 
known  about  the  statistical  standard   deviation  of  Y- ,   each 
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W   can  be   equal  to  one  and   the  S  can  be   determined  by  the 
trial  and  error  method. 

1.   Interpolation   and    Approximation   Using    B-spline 
Function 

There  is  a  difference  between  the  interpolation  and 
the  approximation.  In  interpolation,  the  number  of  knot 
positions  required  are  equal  to  the  number  of  sampling 
points  and  the  value  of  S  in  eq(4.7)  is  small.  In  this  case, 
the  computation  time  required  will  increase  tremendously. 
Whereas,  in  approximation,  it  is  not  of  great  concern  that 
the  approximated  value  at  a  position  be  the  same  value  of 
the  sampling  point,  and  most  of  the  information  still  remain 
in  the  original  curve.  Thus,  decreasing  the  value  of  S  will 
result  in  the  large  number  of  the  knot  positions  and  the 
final  curve  obtained  will  be  similar  to  the  original  curve. 

The  justification  for  using  approximation  rather 
than  interpolation  approach  is  that  though  the  resulting 
curve  may  not  be  the  best  fitting  curve,  it  is  smooth  and 
close  enough  to  the  original.  The  approximation  spline  func- 
tion has  less  number  of  knots  than  the  number  of  samples, 
which  reduces  the  total  processing  time.  In  approximation 
approach,  when  the  appropriate  choice  of  S  value  is  made,  it 
will,  in  turn,  generate  sufficient  number  of  knot  positions 
required  to  provide  a  close  approximation  to  the  original 
curve.  The  ratio  of  sampling  point  to  spline  coefficient  is 
10  to  1  as  shown  in  Figure  4.1.  In  Figure  4.1,  a  guided 
missile  cruiser  ship  at  a  distance  of  45000  feet,  with  290 
sampling  points  and  the  associated  B-spline  knots  are  shown. 
The  original  samples  are  plotted  in  solid  curve.  The  knot 
position  are  show  by  small  circles  in  Figure  4.1.  After 
B-spline  approximation,  the  number  of  the  knots  are  reduced 
to  33.  Hence,  this  technique  is  essentially  a  kind  of  data 
compression  technique. 
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Figure  4.1 


Plot  of  the  Original  Data  and  the  Approximate 
with  Free  Knot. 


TO   DETERMINE 
COEFFICIENTS 


THE   KNOT   POSITION   AND 


THE   B-SPLINE 


The  approximate  smoothing  factor  in  eq(4.7)  is  to  be 
calculated  before  using  the  subroutine  "PARAM"  [Appendix  C]. 
In  the  subroutine,  there  is  a  check  on  S  factor  every  time 
it  is  run.  If  the  S  input  values  do  not  satisfy  or  meet  the 
criteria,  the  program  will  return  with  a  error  code  IER. 

There  is  a  parameter  NEST  which  is  set  to  a  constant 
value.  This  value  determine  the  dimension  of  the  knot  posi- 
tions which  relates  to  the  array  T(NEST),  Cx(L.NEST)  and 
Cy(L.NEST).  The  NEST  is  an  overestimate  of  the  dimension  of 
the  arrays  set  by  the  user.  The  limitation  on  the  value  of 
NEST  for  subroutine  "PARAM"  are  as  follows  [Appendix  C] 

1.  2k+2    NEST    M+k+1 

2.  Typically,  value  of  NEST   M/2 

where  M  is   the  number  of  the  total  sampling   points  and  k=3 
is  the  highest  order  in  the  B-spline  function. 

The  subroutine  "PARAM"  produces  several  outputs  as,  N 
the  total  number  of  knot  positions,  T(N)  the  array  of  the 
value  of  knot  positions  in  Z  parameter,  Cx(N)  and  Cy(N), 
B-spline  coefficient  of  X  functions  and  Y  functions  for  knot 
positions  defined  in  array  T(N) . 
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1.   Limitation  on  B- spline  Coefficient  Determination 

If  the  value  of  NEST  is  too  small,  the  user  will 
receive  error  code  IER  and  the  number  of  knot  positions  will 
appear  to  be  dense  in  the  first  part  of  the  curve.  While 
sparse  in  the  other  part.  The  example  of  an  incorrect  selec- 
tion -of  NEST  is  shown  in  Figure  4.2. 


Figure  4.2    Knot  Selection  if  the  NEST  Parameter  is  too  Small 

Another  problem  which  causes  difficulty  in  program- 
ming is  that  the  main  program  is  in  PASCAL  while  the  subrou- 
tine is  in  FORTRAN.  As  for  the  PARAM  program  in  FORTRAN,  the 
array  index  value  starts  from  1,  while  ,  for  the  user  PASCAL 
program  it  starts  from  0.  Therefore,  in  linking  the  main 
program  to  the  subroutine,  one  has  to  keep  in  mind  of  the 
difference.  In  addition,  the  FORTRAN  programming  logic 
structure  is  so  complicated  that,  when  an  error  occurs,  it 
is  very  difficult  to  debug  and  locate  the  error. 

The  values  of  Cx  and  Cy  depend  upon  uneven  knot 
positions,  and  they  contribute  controlling  effect  to  the 
reconstructed  curve  which  would  be  both  smooth  and  close  to 
the  original  curve.  In  running  the  B-spline  approximation 
program  for  the  first  time,  the  values  of  Cx  and  Cy  of  the 
last  4  knots  at  the  right  end  are  zeros.  However,  in  running 
it  again,  increase  the  value  of  S  did  not  yield  zero  values 
of  Cx  and  Cy  at  those  points.  This,  nevertheless,  has  no 
effect  on  the  reconstruction  of  the  curve. 
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Plots  of  the  reconstructed  profile  of  three  ships  of 
the  CGN  class  using  different  value  of  S,  result  in 
different  number  of  knot  positions  and  are  shown  for  compar- 
ison in  Figure  4.3,  4.4  and  Figure  4.5.  The  original 
samples  are  plotted  in  solid  curve,  the  reconstructed  curve 
is  plotted  in  dashed  line  in  Figure  4.3. 
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Figure  4.3    A  CGN  at  Range  of  45000  feet 
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Figure  4.4    A  CGN  at  Range  of  55000  feet 
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2.   Selection  of  Smoothing  Factor  Value(S) 

It  is  found  by  trials  and  errors  using  the  computer 
program  PARAM  that  in  order  to  retain  maximum  information  of 
the  profile  curve,  the  number  of  knot  positions  required 
should  be  in  the  range  of  25  to  35.  The  number  of  knot 
positions  depends  upon  the  value  of  S  which  must  be  set 
accordingly.  The  appropriate  choice  of  S  to  satisfy  the 
condition  stated  above,  is  important.  For  the  class  of  a 
guided  missile  cruiser (CGN) ,  three  ship  images  at  3 
different  ranges  were  selected.  Then,  run  the  appropriate 
program  for  various  S  factors  to  see  how  the  number  of 
knot(N)  will  vary.  The  results  are  shown  in  Figure  4.6. 
Plots  of  N  vs  S  in  Figure  4.7  through  Figure  4.13  show  that, 
in  most  cases,  the  value  of  N  decrease  quite  rapidly  when 
the  value  of  S  is  in  the  range  of  0  to  100,  and  gradually 
for  S  factor  in  the  range  of  100  to  200,  thereafter,  the 
value  of  N  decreases  very  little.  Obviously,  the  curve  seems 
to  decay  exponentially.  Furthermore,  for  some  classes  of 
ships  changes  are  more  pronounced  than  the  others  which  is 
probably  due  to  the  actual  number  of  knots  present  in  the 
profile.  For  guided  missile  cruiser  ship(CGN)  with  2  lumps, 
the  number  of  knot  positions  required  can  be  33  as  shown  in 
Figure  4.6  We  select  the  factor  S  to  be  about  100. 

The  selection  of  the  factor  S  depends  upon  the 
number  of  the  original  sampling  points.  If  the  factor  S  is 
small,  large  number  of  knots  are  needed.  When  the  factor  S 
is  large,  small  number  of  knots  are  needed.  When  the  number 
of  knot  and  the  Cx  and  Cy  coefficients  are  small,  the 
B-spline  coefficients  Cx  and  Cy  obtained  can  not  be  used  to 
reconstruct  the  curve  close  to  the  original  profile. 
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3.   The  Output  Cx,Cy_  and  Knots  Profiles 

From  the  previous  study  the  factor  S  is  selected  for 
each  of  the  class  of  ship  as  follows 

1.  Destroyer(DD) ,  S  =  20.0 

2.  Container,  S  =  100.0 

3.  Freighter,  S  =  100.0 

4.  Replenishment  oiler,  S  =  20.0 

5.  Tank  landing  ship(LST),  S  =  25.0 

6.  Frigate(FF),  S  =  100.0 

7.  Guided  missile  cruiser ( CGN) ,  S  =  125.0 

8.  Guided  missile  destroyer (DDG) ,  S  =  125.0 

The  plot  of  the  B-spline  coef f icients , Cx  and  Cy, and 
X,Y  at  the  positions  of  the  sampling  points  vs  the  Z  param- 
eter are  shown  in  Figure  4.14  through  Figure  4.21. 
Observation  and  comparisons  of  the  curves  show  that  the 
values  of  Cx  and  Cy  exhibit  changes  similar  to  that  of  X  and 
Y  except  that  the  variation  of  values  leads  that  of  the  X 
and  Y.  This  is  due  to  the  fact  that  Cx  and  Cy  have  to  act 
as  the  controlling  factor  for  the  reconstructed  B-spline 
curve  to  get  the  result  close  to  the  original  curve. 

Examination  of  the  plots  of  the  results  of  X  and  Y 
show  that  when  X  is  increasing  monotonically ,  Y  is  almost 
constant;  but  when  X  is  almost  constant,  Y  is  increases  or 
decreases.  This  behavior  relating  to  profile  reconstruction 
may  be  explained  as  follows.  For  a  ship  profile  when  X  is 
increasing  and  Y  not  increase  too  much,  this  may  be  inter- 
preted as  an  almost  leveled  profile.  When  X  is  almost 
constant,  and  Y  may  be  increasing  or  decreasing,  it  may  be 
interpreted  as  the  beginning  or  the  ending  of  the  lump. 
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D.   SHIP  DESCRIPTION 

The  shape  of  the  ships  depend  upon  the  shapes  and  posi- 
tions of  the  lumps.  Characteristics  of  the  lumps  for 
different  type  of  ships  are  as  follows: 

1.  Frigate  -  The  beginning  of  the  lump  is  at  1/6  of  the 
ship  length  to  the  left  side  of  the  midships  and 
the  peak  of  the  lump  is  at  the  mast  which  is  located 
at  the  midships.  The  termination  of  the  lump 
is  approximate  1/3  of  the  ship  length  from  the 
midships  to  right  side.  In  addition,  the  average 
height  of  the  lump  is  a  little  higher  than  the 
level  between  the  bow  and  the  stern,  and  its  size 
is  1/2  of   the  ship  length  as  shown  in  Figure  4.22. 


GASClA  Ff  (with  lampi) 


USA  18) 


Figure  4.22    Frigate. 

Container  -  The  beginning  of  the  lump  is  at  1/3  of 
the  ship  length  to  the  right  side  of  the  midships 
while  the  lump  is  high  and  terminate  at  the  stern. 
The  lump  appears  to  be  in  a  rectangular  shape  with 
small  crane. 

Tank  landing  ship(LST)  -  The  beginning  of  the  lump  is 
at  1/4  of  the  ship  length  from  the  midships  to  the 
left  side;  the  height  of  the  lump  is  higher  than  the 
level  between  the  bow  and  the  stern  by  a  small 
margin,  while  its  highest  point  is  located  at 
approximate  1/6  of  the  ship  length   from  the  midships 
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to   the  left   side   with   small   difference   from  the 
average  high  as  shown  in  Figure  4.23. 


NEWPORT.  -NfWPORT-   Cl«£i  LST 


USA  :20l 


Figure  4.23    Tank  Landing  Ship(LST) 


6. 


Guided  missile  destroyer  -  the  beginning  of  the  lump 
starts  at  approximately  1/4  of  the  ship  length  from 
the  midships  to  the  left  side  with  its  highest  point 
at  the  mast.  After  this  the  slope  will  decline  in  a 
very  rapid  fashion  following  a  noticable  deck 
distance,  then  the  beginning  of  the  second  lump 
occurs  due  to  the  redome  presence.  Therefore,  the 
lump  will  be  narrow  with  great  height,  and  will 
terminate  at  approximatly  1/6  of  the  ship  length  from 
the  midships  to  the  right  side.  Furthermore,  there 
is  a  little  lump  near  the  stern,  this 
distinguish  destroyer  of  the  same  size  as  shown  in 
Figure  4.24. 

Destroyer-lump  characteristic  will  be  similar  to  that 
of  guided  missile  destroyer  except  for  the  small  lump 
as  in  Figure  4.25. 

Guided  missile  cruiser  -  The  beginning  of  the  first 
lump  is  at  1/12  of  the  ship  length  from  the  midships 
to  the  left  side  with  highest  point  at  the  mast.  The 
lump  size  is  large  both  in  length  and  height  and  its 
height  decreases  to  the  point  which  is  a  little  above 
the  level  between  the  bow  and  the  stern.  Therefore, 
the  second  lump   begins  with  almost   the  same  size  as 
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Figure  4.24    Guided  Missile  Destroyer (DDG) 


Figure  4.25    Destroyer. 

the  first  one.  The  second  lump  ends  at  approximately 
1/4  of  the  ship  length  from  the  midships  to  the  right 
side.  Furthermore,  the  distance  between  the  peak  of 
both  lumps  is  less  than  that  of  the  Guided  missile 
destroyer  shown  .in  Figure  4.26. 


BAINBRIDGE    "BAIN8RI0GE"  Class  CGN 


USA  111 


Figure  4.26    Guided  Missile  Cruiser(CGN) 
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7.  Replenishment  oiler  (AOR)  -  There  are  2  little  lumps 
with  the  first  one  starting  at  4/5  of  the  ship  length 
from  the  midships  to  the  left  side  exhibiting  rapid 
decrease  in  height  to  the  level  between  the  bow  and 
the  stern.  The  second  lump  starts  at  1/4  of  the  ship 
length  from  the  midships  to  the  right  side  with  slow 
decline   to   the   point   near   the  stern. 

8.  Freighter-with  3  noticeable  lumps,  the  first  two 
represent  lifting  crane  with  the  first  high  lump 
approximately  1/6  of  the  ship  length  from  the 
midships  to  the  left  side  but  narrow;  meanwhile,  the 
second  lump  is  located  at  approximately  the 
midships  with  the  same  size  as  the  first  one.  The 
beginning  of  the  third  is  at  approximately  1/3  of 
the  ship  length  from  the  midships  to  the  right 
side  with  large  size.  It  is  higher  than  the  first 
two  and   with  gradual  decline  to  the  stern. 

From  the  lump  characteristic  of  different  type  of  ships, 
we  could  distinguish  the  type  of  a  ship  based  on  the  rotated 
superstructure . 

E.   SHIP  CLASSIFICATION 

An  important  aspect  in  categorizing  types  of  ships  is  to 
recognize  the  lumps  above  the  main  deck.  Therefore,  it  is 
useful  to  plot  the  positions  of  X  and  Y  vs  the  Z  parameter 
where  the  Z  parameter  can  be  determined  from  eq(4.5)  where 
Z(I)  is  an  array  (1..M)  as  shown  in  Figure  4.27  where  X  is 
ploted  in  solid  line  and  Y  is  ploted  in  dash  line.  Then, 
plot  the  B-spline  coefficients  Cx  and  Cy  vs  T(N);  the  knots 
position  in  Z  where  N  is  the  total  number  of  knots  as  shown 
in  Figure  4.28.  The  Cx  is  ploted  in  solid  line.  The  Cy  is 
ploted  in  dash  line.  The  Values  of  Cy  will  be  close  to  the 
curve  of  Y   while  the  values  of   Cx  for  the  same   knot  posi- 
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tions  as  Cy  is  exhibiting  monotonic  increasing  trend.  When 
the  difference  in  Cx  is  decreasing  or  zero,  the  value  of  Cy 
will  have  a  pronounced  change  where  the  beginning  or  the 
ending  points  of  a  lump  can  be  detected.  The  value  of  knot 
position(T)  at  those  points  related  to  the  sizes  of  the 
Cx,Cy  can  be  determined.  Finally,  with  the  values  of  Cx  and 
Cy  at  those  points  known,  the  area  of  the  lumps  can  also  be 
determined. 

Hence,  from  the  ship's  characteristics  and  information 
derived  from  the  above  procedure,  classification  of  ships 
can  be  made  by  considering 

First,  the  number  of  lumps  detected,  1,  2,  or  3 . 

1.  The  number  of  lump=l:  Frigate,  Tank  landing  ship 

2.  The  number  of  lump=2 :  Destroyer,  Guided  missile 
cruiser,  Guide  missile  Destroyer,  and  Replenish- 
ment oiler (AOR) 

3.  The  number  of  lump=3 :  Freighter  and  Container 
Second,   the  position  of  a  lump  relative  the  midships  is 

measured.  This  quantity  is  scaled  by  the  total  length  of 
the  ship.  This  scaled  quantity  will  be  invariant  with 
respect  to  the  different  ship  sizes  at  different  ranges 

Third,  the  area  of  the  lump  is  normalized  to  the  ship 
length  squared. 

1 .   Lump  Detection 

As  shown  in  the  plot  of  X,  Y,  Cx,  and  Cy  vs  Z  param- 
eter in  Figure  4.15,  when  the  difference  (ACx)  between 
successive  value  of  Cy  varies  from  increasing  to  decreasing, 
and  then,  to  increasing  again,  Cy  exhibits  noticable  varia- 
tion. From  observation  of  Cx  values,  it  is  seen  that  the 
difference  (ACx)  always  has  variation  in  the  same  sequence 
as  stated  above.  Therefore,  only  Cy  values  are  taken  into 
consideration  in  the  program  that  detects  lump. 
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There  are  two  distinct  procedures  used  in  dealing 
with  two  different  kinds  of  lumps  big  or  small.  Thus,  it  is 
necessary  to  distinguish  in  the  first  place  whether  the  lump 
is  big  or  small.  For  the  big  lump,  the  differences  of  Cy  is 
positive  for  all  initial  four  knots.  It  continues  to  the 
peak  and  decreases  toward  the  ending  of  the  lump.  For  the 
small  lump,  the  difference  of  Cy  is  positive  for  the  first 
two  knots  and  stay  constant  or  positive  for  the  third  knot, 
but  the  difference  of  Cy  for  the  forth  knot  is  negative, 
thereafter,  the  program  proceeds  to  establish  the  following 
values  for  each  lump  detected: 

First,  the  knot  positions  (Z  value)  at 

1 .  the  beginning  of  the  lump 

2 .  the  ending  of  the  lump 
Second,  the  Knot  number  for 

1 .  the  beginning  of  the  lump 

2.  the  ending  of  the  lump 
Third, Number  of  lumps  detected. 
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a.   Procedure  to  Detect  the  Lump 

1.  Check  a  big  or  small  lump  by  testing  the  conditions 
of  three  varying  values  of  Cy  increments  (A  Cy)  in 
sequence.  If  they  are  positive  it  is  a  big  lump  and 
set  the  flag- lump  to  1,  otherwise  it  is  a  small  lump 
and  set  the  flag- lump  to  zero;  then  go  to  2 

2.  Check  the  present  knot  position  to  see  where  its 
Z  value  equal  to  zero  or  to  maximum  Z  value,  if  it  is 
Z  maximum  then  stop,  it  not  go  to  3 

3.  If  the  process  begins  at  the  first  knot  position,  set 
the  begin  and  the  flag-end  to  zero;  check  the  status 
of  the  flag- lump  for  1  (big  lump)  or  0  (small  lump), 
if  it  is  a  big  lump,  then  go  to  4.  If  it  is  a  small 
lump,  then  go  to  7. 

4.  Check  the  status  for  the  begining  or  ending  of  the 
lump.  If  the  flag-begin  is  1,  it  represents  that  the 
beginning  of  a  lump  is  found,  then  go  to  5.  Otherwise 
it's  not  found,  then  go  to  6. 

5.  Find  the  ending  of  the  big  lump  by  testing  the 
conditions  of  3  values  of  Cy  increments  in  sequence. 
The  first  2  should  be  negative  and  the  third  should 
be  constant  or  positive.  If  the  condition  are 
satisfied,  store  the  Z  value  of  the  position  of  the 
third  knot,  and  the  flag-begin  to  zero;  then  go  to 
10. 

6.  Find  the  beginning  of  the  big  lump  by  testing  the 
conditions  of  3  different  values  of  Cy  increment 
( A  Cy)  in  sequence.  The  first  Cy  should  be  negative 
or  constant,  the  second  should  be  positive,  and  the 
third  should  be  positive.  If  the  conditions  are 
satisfied,  store  the  Z  value  of  the  position  of  the 
second  knot,  and  set  flag-begin  to  1;  then  go  to  10. 

7.  Check  the  status  for  the  beginning  or  ending  of  the 
lump.   If   the  flag-begin  is   1,   it   represents   the 
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beginning  of  a  lump  is  found,  then  go  to  8.  Otherwise 
it's  not  found,  then  go  to  9. 

8.  Find  the  ending  of  the  small  lump  by  testing  the 
conditions  of  3  values  of  Cy  increment  (A  Cy)  in 
sequence.  The  first  one  should  be  negative  or 
constant.  If  the  conditions  are  satisfied,  store  the 
Z  value  of  the  position  of  the  second  knot,  and  set 
the  flag-begin  to  zero;  then  go  to  10. 

9.  Find  the  beginning  of  the  small  lump  by  testing  the 
conditions  of  3  values  of  Cy  increment  (A  Cy)  in 
sequence.  The  first  should  be  constant  or  negative, 
the  second  should  be  positive,  the  third  should  be 
constant.  If  the  conditions  are  satisfied,  store  the 
position  of  the  second  knot,  and  set  the  flag-begin 
to  1;  then  go  to  10. 

10.   Move  to  the  next  knot,  then  go  to  2. 
This  procedure   is  shown   in  Figure  4.29   and  the   detail  of 
each  procedure  is  shown  in  Appendix  D. 
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2 .   To  Determine  the  Area  Under  a  Lump 

The  area  under  the  lump  can  be  determined  by 


AAREA  =  ACxACy  +  CyACy 
2 


(4. 10) 


Suppose  there  is  a  lump  as  shown  in  Figure  4.30_. 
The  area  increment  can  be  calculated  by  using  eq(4.10).  Then 
the  area  under  BC(which  is  negative)  is  added  to  the  area 
under  AB . 
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y         ^ 
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p 

/3 

/+ 

».Pv 
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•>      Q             E 

"■  oX 

Figure  4.30    First  Procedure  to  Determine  the  Area. 

Next,  area  under  CD  is  calculated,  which  is  positive 
and  add  this  to  the  last  resulting  area.  Obviously,  the  sum 
would  represent  the  total  area  of  the  lump  as  shown  in 
Figure  4.31. 
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Figure  4.31    Step  by  Step  Procedure  to  Determine  the  Area. 

F.   CONSTRUCTION  OF  THE  DECISION  TREE 

The  characteristics  of  ships  used  for  classification 
are,  the  number  of  lumps,  the  area  of  lumps,  the  knot  posi- 
tions (Z  value)  and  where  the  beginning  of  the  lump  and  lump 
maxima  relative  to  the  midships  are  located. 

In  view  of  the  above  characteristic,  it  is  necessary  to 
find  the  relationships  among  them  to  make  classification 
possible.  Therefore,  for  each  of  the  eight  different  class 
of  ships,  the  decision  tree  is  constructed  which  is  based 
upon  relationships  observed  in  the  plots  of  the  knot  posi- 
tion (Z  values)  for  the  beginning  of  the  lumps  normalized  by 
the  total  ship  length  (Z  value)  VS.  the  number  of  lumps;  the 
area   of   the  lumps   normalized   by   the  total   ship   length 
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(Z  value)  squared  vs  the  number  of  lumps;  and  the  Z  value  of 
the  lump  maxima,  as  shown  in  Figure  4.32  through 
Figure  4.39. 

Thus,  according  to  the  number  of  lumps  presented  in  the 
profile  used  as  the  first  criteria,  ships  can  be  divided 
into  3  distinct  groups.  Then,  classification,  can  be  made 
by  comparing  further  characteristic  as  shown  in  Table  III, 
and  the  complete  decision  tree  constructed  from  Table  III  is 
shown  in  Figure  4.40. 
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Figure  4.32    Plot  the  Beginning,  Area,  and  Peak  of  a  FF 
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Figure  4.33    Plot  the  Beginning,  Area,  and  Peak  of  a  LST 
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Figure  4.34    Plot  the  Beginning,  Area,  and  Peak  of  a  DD 
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Figure  4.35    Plot  the  Beginning,  Area,  and  Peak  of  a  DDG 
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Figure  4.36    Plot  the  Beginning,  Area,  and  Peak  of  a  CGN 
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Figure  4.37    Plot  the  Beginning,  Area,  and  Peak  of  a  AOR. 
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Figure  4.38    Plot  the  Beginning,  Area,  and  Peak  of  a  Freighter 
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TABLE  III 
Comparison    of    Different   Types    of    Ships 


|No.    OF    LUMP 

EEC  INN  INC   OF    LUMP 
FROM   MIDSHIP 

PEAK    IN 

THE    LUMP 

AREA   WIDER    THE    LUMP  | 

RESULT 

1 

-0.4«BECIN«-0.3 

-    0.2s:FEAK*0.1 

0.0C5<AREA<0. 007         | 

FF 

-0.2i5ECINtf0.0S 

-0.  1*PEAK*0 

O.C01*AREA*0.005         | 

LST 

-0.3<lst    EECIN*-0.2 
0     <•  2nd    BECIN*0.  1 

-0.2*lst 
0.05<2nd 

PEAK*-0.  1 
PEAK^O. 15 

0.002*lst 
.002S*2nd 

AREA*0.00S| 
AREA*,.  0055| 

CCN 

2 

-0.2*lst    BEGIN*-0.05 
0.  1*  2nd    BEGIN*0.2 

-0.1*lst 
0.2*;  2nd 

PEAK<0 
PEAK*0.3 

0.002*lst 
0.001*2nd 

AREAS  0.0031 
AREA*j.CC25l 

DDG 

-0. 4*lst    BECINtfO. 1 
-0. 5*2nd    BECIN<0. 1 

-0.4*lst 

0.  l<:2nd 

FEAK*0 
PEAKtf0.25 

.0005«lst 
0. 004<2nd 

AREA*0.002| 
AREAS0.007I 

AOR 

-C.5*lst    BECIN*-0.3 
-.15<2nd    EECIN*0.0S 

-Q.Zclst. 
0    <  2nd 

PEAK* 0.05 
PEAKtf.0.1 

0.002*lst 
0.001*2nd 

AREA*0.008| 
AREA40.007] 

DD 

|            3 

-0.5*lst    BEGINc-0.4 

-0.2<,2nd    BECIN<-0.  1 

0. K3rd    BECIN*0.2 

-0.5*lat 

-0.1<2nd 

0.3Ord 

PEAK* -0.4 
PEAK*    0 
PEAK* 0.4 

0.001*lst 
0. 001* 2nd 
0. 003<3rd 

AREA*0.002| 
AREA*0.002|c 
AREA<0.  004| 

REIGHTER 

-0.3*lst    BECIN*-0.2 
-0. l*2nd    BECIN<  0 
0.  15s3rd    BECIN<0.25 

-.15* 1st 

-  . 05*2nd 
0.2*3rd 

PEAK*- .05 
PEAK*D.OS 
PEAK*;0.3 

0     41st 

0    -£2nd 

.003543rd 

AREA*0.001| 

AR£A*0.  00 1|  CONTAINER 

AREA*.0045| 
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G.   SUMMARY 

The  decision  tree  which  is  constructed  from  Table  III 
does  not  provide  100-percent  correct  classification  for  all 
types  of  ship  images.  Furthermore,  the  presence  of  noise  in 
images  may  cause  complicatations  in  classifying  the  ships. 
For  example,  with  excessive  noise  in  the  original  image, 
Sobel  operator  for  edge  enhancement  in  the  preprocessing 
process,  still  yield  result  with  residual  noise  present  in 
Figure  4.41.  These  residual  noise  is  undesirable  since  it 
causes  failure  to  extract  profiles  which  retain  necessary 
informations  from  the  original  images.  The  subsequent  clas- 
sification of  profile  by  the  Fourier  Coefficient  method  and 
the  B-spline  Coefficient  method  becomes  difficult. 


Figure  4.41    Noisy  Image. 

As  discussed  before  in  order  to  reduce  the  noise,  an 
appropriate  threshold  gray  value  for  the  preprocessed  image 
is  set.  This  results  in  a  silhouette  image.  However,  if  the 
value  is  too  high  the  profile  becomes  broken  which  prevent 
successful  operation  in  the  closing  process  discussed  in 
chapter  2,  as  shown  in  Figure  4.42.  On  the  other  hand 
decreasing  the  threshold  value  results  in  erroneous  profile, 
as  shown  in  Figure  4.43.  In  some  cases,  when  the  closing 
process  is  used,  certain  vital  information  is  lost.  This 
effect  can   be  seen  in   comparing  the  image   in  Figure  <*.44, 


100 


and  Figure  4.45.  Therefore,  loss  of  informations  due  to  the 
attempt  of  eliminating  noise  and  the  inability  in  the 
preprocessing  process  to  cope  with  the  residual  noise,  yield 
erroneous  profile.  Consequently,  failure  occurs  in  the 
succeeding  classification  steps. 


Figure  4.42    High  Threshold 


Figure  4.43    Low  Threshold 
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Figure  4.44    Before  Closing 


Figure  4.45    After  Closing 
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V.  CONCLUSION 

The  shape  of  the  superstructure  of  a  ship  is  the  most 
important  feature  used  in  the  classification  algorithms.  To 
extract  the  edge  profile,  a  Sobel  operator  is  employed.  The 
limitation  of  a  Sobel  operator  is  that  some  small  details  of 
the  edge  is  lost.  For  example,  the  image  of  a  DD  at  a  range 
of  79000  feet  after  applying  the  Sobel  operator  shows  the 
superstructure  and  the  radar.  But  the  edge  of  a  small  mast 
disappeared  as  shown  in  Figure  5.1.  Furthermore,  if  the 
threshold  value  is  set  too  big,  the  edge  image  of  the  super- 
structure profile  becomes  broken.  The  top  superstructure 
profile  is  obtained  by  setting  the  gray  value  under  the 
slope  between  the  bow  and  the  stern  to  zero.  We  need  to 
apply  a  contour  tracking  process  to  refine  the  superstruc- 
ture profile.  For  some  images,  the  connection  of  the  broken 
profile  pieces  may  be  achieved  in  a  Closing  operation  as 
discussed  in  Chapter  2.  The  disadvantage  of  the  Closing 
operation  is  that  some  small  details  of  the  profile  may 
disappear.  The  superstructure  profile  of  a  freighter  at  a 
range  of  53000  feet  is  shown  in  Figure  5.2.  After  the 
Closing  process  the  detail  of  the  cranes  disappeared  as 
shown  in  Figure  5.3. 

The  ship  classification  can  be  achieved  by  either  the 
Fourier  Coefficient  method  or  the  B-spline  Coefficient 
method.  In  using  these  coefficients  to  classify  ships,  it 
is  found  that  only  the  initial  coefficients  that  lie  between 
the  0th  to  the  20th,  are  relevant,  while  the  rest  are  not. 
Inspection  of  the  comparison  curves  of  the  same  class  of 
ships  shows  that,  similarities  in  patterns  exists  up  to  the 
20-th  point.  Beyond  that,  diversities  in  shape  are  so  great 
that  inclusion  of  those  additional  points  for  classification 
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will  be   of  no  use. 
Figure  5.5. 


Examples   are  shown  in   Figure  5.4  and 


Figure  5.1    Edge  Image  of  a  DD  at  a  Range  of  79000  feet 


Figure  5.2    Before  Closing  Process 


Figure  5.3    After  Closing  Process 
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In  the  B-spline  coefficient  method  the  technique  of 
uneven  knot  selection  has  been  employed.  With  the  original 
sampling  points  on  the  ship  profile  as  input,  a  smaller  set 
of  approximation  samples  are  collected.  These  can  be  used  to 
reconstruct  the  ship  profile  with  sufficient  information 
retained  from  the  original  image.  The  execution  time  as 
approached  to  that  of  handling  the  original  data  directly 
has  been  greatly  reduced.  The  reduction  of  the  number  of 
sample  points  by  almost  a  factor  of  10  is  common.  For  a  CGN 
ship  a  set  of  original  sampling  points  of  290,  has  been 
reduced  to  36. 

Comparing  the  two  methods,  Fourier  Coefficient  and 
B-spline  Coefficient,  the  former  method,  in  some  cases  is 
not  effective  to  establish  satisfactory  classification  of 
ships.  This  is  due  to  difficulties  in  matching  similarities 
of  the  shape  of  the  coefficient  curves.  The  latter, 
however,  surpasses  the  former  in  that  it  is  able  to  classify 
more  ships  accurately  using  computer  programs.  It  is 
possible  to  improve  the  reliability  of  those  two  methods  by 
reducing  noise  in  the  data  collection  process  and  the 
preprocessing  process. 
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APPENDIX  A 
THE  PROGRAM  TO  OBTAIN  THE  SUPERSTRUCTURE 
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Source    Listing 

program    cut(input,outputunfile,outfil9)  ; 

type 

byte    =    0..255I 

lmagarouil     =    packed     array    CO.. 257]    of     byte; 

imagarou/2     =     pack°d     array     CO.  .2552     of     lnteoer; 

imag?roui3    -    packed     array 

roul        -     oacked     array     CO.. 


6-Dec-198<.  17:18:40 
6-Dec-198*  17:18:31 


VAX- 
_D.RA 


CO. .2553  of  byte: 
,255]  of  byte". 


ro«2   =  oacked  array  CO. .255]  of  integer; 


sooel  :  array  CO. .63]  of  inagerow2: 

it  J   :   byte: 

f  :  array  CO. .65]  of  imagerowi: 

outfile  :  file  of  imager  oiu  3 ; 

d« t dy  ,  range . or ight » n , m     lintegeri 

infile  :  file  of  rovi; 

slope  :  real; 

image  larray    CO. .633  of  roml'. 

x  i  x 1  i  x2 i y 1 i yZ     '     integer; 

NUM1 ,NUM2 ,NUM3,NUM4, max ,min, maul  :  int?9er; 

cal  larray     CO. .633  of  roul; 

thes  :  integer: 

NAME  :  PACKED     ARRAY  CI. .203  OF  CHAR: 

BEGIN 

JRITELNC  'INPUT  SHIP   FILENAME  OUT  CUT. DAT'): 

READLNCNAME) ; 

MR  ITELNC  'THRESHOLD*): 

READCTHES) : 

WRITELNC'NUM  TC  CUT  FROM  LEFT*): 

REAO(NUMl): 

URITELNCNUM2  TO  CUT  FROM  RIGHT*): 

READCNUM2) : 

WRITELNC 'NUM3  TO  CUT  FROM  TOP'): 

REA0CNUM3) ; 

WRITELN(  'NUM<»  TO  CUT  FROM  BUTTOH'); 

REA0CNUM4) : 

open  ( infile f NAME , hi  story  :  =  old, 
ac cess_me t hod  : ■ sequen t ial > 
r ec or d _leng t h  : = 25 6 t r ec ord _ t y pe  :=fixed): 

open  ( out f i 1 e , *C UT . dat * , hi s t or y  : =new » record _ 1 enqt h 
record_type  :=fixed): 

reset(infile)  ; 

rewrite  (outfile); 

i:  =  o; 

»hile  not  eof  (infile)  do 

oegin 

read  ( in f l le , imageCi 3)  I 
for  j:=0  to  255  do 

fCi*ltj*13  :=  imageCitj]* 

i:=i*l : 

end  ; 

(^compute  sooel*) 


=  256, 
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Source  Lasting 


6-0ec-198<. 
6-Dec-1984 


17:18:40 
17:18:31 


VfiX- 
_DTR< 


FOR  j:=0  TO  7  00 

cC64,J  +  2>.Sj:  =  FC63,J+2483: 
for  i  :=  0  to  63  ao  begin 

fti*l,03  :=  fCi-H.13: 

fCi+1.2573  :=f Li*l,2563; 
«nd : 
for  j:=  0  to  255  do  begin 

fco.j-na  :=  fcitj+13: 

fC65,j*13  :=  fC64fj*13; 
end  ; 

fC3.03:=fCl, 13: 
fC0,2573:=fCl,2563: 
f C65,03:=f C64.13; 
f C65.2573:=fC64,2563: 

(*set  initial  ia<i  min  *) 

■  si  :=o ; 
max i :  =  o ; 
min:  =  25S  ; 

for  i:=  0  to  63  do  begin 
for  j:=  0  to  255  do  begin 

dxisfCit j3>2*fCi*l. jD*fCi  +  2. j3— fCi» j  +  21 

-2*fCi*l,j«-2  3-fCi  +  2tj*23; 
dy:  =  fCi*2,j3-»2*fCi-*2,j-U3-*fCi*2,j*2  3-fCi,j3 

-2*f Li, j*ll-i Cii j+235 
sobelCitj3:=rojnd(CdK**2*dy**2)**0.5>; 
if  max  <  sobelCitjJ  than 

max  :=  sobelCi.j}; 
if  min  >  sodeKii  j]  then 
min  :=  sobelC i t j 3 : 
end  ; 
end  ; 
range  :=  max-min; 

(*  rescale*) 

for  i:=  0  to  63  do  begin 
for  j:=  0  to  255  do  begin 

calCiij3  :=  roundC ( C sobelCi > j3-min)*255 )/r ange ) 
if  (j<=NUMl)  or  ( J>=255-NUM2)  then 

calCi, j3  : =0  : 
if  (i<=NUM3)  or  (i>=63-NUKO  then 
calCi,j3  :=  0: 

(*  cut  thesholdO 

IF  THESO0  THEN  BEGIN 
if  calCi,J3  <-th9S 
then 

calCi,j3  :=  0 
else 

calLi,j3  :=  255: 

eno; 

end: 
end: 


(*find  point  at  raw  and  aft*) 
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Source  Listing 


6-0ec-198'i 
6-£jec-1964 


17:i3:<.0 
17  :1B  :31 


VAX- 

_DRA 


m  :=  o : 

n  :  =  o: 

bright:=  255: 

for  j  :=  0  to  150  do  begin 
for  i:=0  to  63  do  begin 

if <calCi i jD=bright  )  and  (n  =  0)  then 

begin 

yl  :=  j: 

xi  :=  i: 

n  : =  n* 1  ; 
end;  (*ifO 

if (calCi, 255- j3=bright)  and  (m=0)  then 
begin 

y2  :=  255-j; 

k2  :=  i: 
n  : =  m* l ; 
end;  (tift) 
end: (*for*) 
end ; ( *f  or*  ) 
■riteln("xl=",xl,*yl=*,yl,*x2="tx2i*y2="»y2) 

C*  cut  line*) 

s  looe  :  =  1  .  0 ; 
if  x 1 s« 2  then 
begin 

s  looe  : =  0 : 

for  i:=xl*l  to  63  do  begin 
for  j:=yl-5  to  y2*5  do 
calti, j3  :=0; 
end ; <  *f  or* ) 
end : (*if *) 
if  slopeOO  then 
begin 

slope  :=(x2-xl)/(y2-yl): 
if  slope  <  0  then  begin 
slope  :  =  abs(slope); 
for  j  :=  yl  to  y2  do  begin 

x  :=  xl  ♦l-roundC  (  j-yl  )<=slope)  : 
for  i  :=x  to  63  do 
calCi  t  j3:  =  0 : 
end  Ktfor*) 
end: (titt) 
if  slope  >0  then  begin 

for  j:=  yl  to  y2  do  begin 

x  :=  x 1 *1 *round( ( j-y 1 )»slope) : 
for  i  :=x  to  63  do 
calCi f j3  :=0: 
end : ( * t  ors) 
end ; (si f «) 
uiriteln("step2"); 
end;(*if«) 

(♦put  outfileO 

for  i:=0  to  63  do  begin 
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S  our c  e  Listing 


6-09C-1 994 
6-Dec-1984 


17:13:40 
17:18:31 


VAX- 


for  j  :=0  to  255  do 

outfiie*Cj3  :=  calCi, 
put  (outfile): 
end : C  *f ort ) 
clossC in  f l le) : 
•  nd. 


j3: 
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APPENDIX  B 
THE  PROGRAM  OF  CONTOUR  FOLLOWING 
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Source  Listing 
PROGRAM  TRACECINPUT.OUTPUT.INFILE.OUTFILE); 

TYPE 

BYTE  =  0..255: 

IMAGEROW1  =  PACKED  ARRAY  CO. .2553  OF  BYTE: 

ROU1  =  PACKED  ARRAY  CO. .2573  OF  BYTE; 

VAR 

R.C,F0,F1,F2,F3.F4,F5.F6,F7,FB   :   BYTE: 

F  :  ARRAY  CO. .653  OF  ROWl: 

A, IMAGE  :  ARRAY  CO. .633  OF  IHAGEROUl; 

INFILE  :  FILE  OF  IHAGEROWl: 

R1,C1,R2,C2.CDIR,I, J, COUNT  :  INTEGER: 

ROW, COL  :  ARRAY  CO.. 5123  OF  INTEGER*. 

OUTFILE  :  FILE  OF  IMAGEROWi: 

NAME  :  PACKED  ARRAY  C1..203  OF  CHAR; 


10-Dec-1984  13:53:33 
6-Dec-198<.  17:33:13 


VAX-11 
QRAO: 


PROCEDURE  STORE: 

BEGIN 

COLCI3:=C-l ; 
ROWCI3:=R-i: 
WRITELN(-COL=*,COLCI3,*ROW=*,RQWCI3): 

count  :=i: 
l:=i*i: 

eno: 


PROCEDURE  CMOVE: 
BEGIN 

fo  :=f:r-i,c3:  fi:=fcr-i,c*13:  f2:=fcr.c*13:  F3:=FCR*itC*i3: 

F<.:=FCR*1,C3:  F5:  =  FCR*1,C-13:  F6:*FCR  ,  C-13  ;  F7  :  =  FCR-1  .  C-13  : 
CASE  COIR  OF 
0:  BEGIN 

IF  (F3=0)  AND  CF2=255)  THEN  BEGIN  C:=C+l;  CDIR:=1:  END 
ELSE  IF  CF2=0)  AND  (Fl=255)  THEN  BEGIN  R:=R-1 

store:  c:=c*i;  cdir:=o:  eno 
else  if  (f1=0)  and  (f0=255)  then  begin  r:=r-1 

cdir:=o:  eno 
else  if  (f0=0)  and  (f7=255j  then  begin  c:=c-1 

store:  r:=r-i;  cdir:=3:  eno 
else  if  cf7=0;  and  (f6=255>  then  begin  c:=c-1 

CDIR:=3:  end 
ELSE  IF  (F6=C)  AND  (F5=255)  THEN  BEGIN  R:«R+i; 

STORE:  C:=C-i;  C0IR:=3;  End 
else  begin  r:=r*i:  c0ir:=2:  eno: 
end: 

l:  BEGIN 

IF  (F5=0)  AND  CF4=255)  THEN  BEGIN  R:=R*i: 

CDIR:=2:  end 
else  if  (f<.=  0)  and  (f3  =  255)  then  begin  c:=c*1; 

store:  r:=r+i:  coir:=i:  end 
else  if  (f3=0)  and  (f2=255)  then  begin  c:=c*1; 

CDIr:=i:  end 

ELSE  IF  (F2=0)  ANC  (Fl=255)  THEN  BEGIN  R:=R-i: 

store:  c:=C*i;  coir:=o;  eno 

ELSE  IF  (Fl=0)  ANC  (FO=255)  THEN  BEGIN  R:=R-1; 
COIR:-0;  ENO 
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Sourc  e  Listing 


10-Dec-1984  13:53:33 
6-Uec-198H  17:38:13 


vax-ll 
_d*rao : 


else  if  (f0=0)  and  <f7=255)  then  begin  c:=c-i: 
store:  r:=r-i;  coir:=o:  end 

else  begin  c:=c-i:  cdir:=3:  end: 
end: 
2:  BEGIN 

if  (f7=0)  and  (f6=255)  then  begin  c:*c-i; 

cdir:=3:  end 
else  if  cf6=0)  and  cf5=255)  then  begin  r:=r*1 

store;  c:=c-i:  cdir:=3:  end 

ELSE  IF  (F5=0>  AND  (F4=255J  THEN  BEGIN  R:=R*1 
CDIR:=2:  ENO 

else  if  cf4=0)  and  (f3=255)  then  begin  c:=c*1 

store;  r:=r*i:  coir:=i:  end 
else  if  (f3=0)  and  (f2=255)  then  begin  c:=c*1 

CDIR:=i;  ENO 
ELSE  IF  (F2=0)  AND  (Fl=255)  THEN  BEGIN  P:=R-1 

store:  c:=c-n:  coir:=i:  end 
else  begin  r:=r-i;  cdir:=0;  eno: 
eno: 

3:  BEGIN 

IF  (F1=0)  ANO  CF0=2553  THEN  BEGIN  R:=R-i; 

cdir:=o;  END 
else  if  cf0=0)  and  (f7=255)  then  begin  c:=:-l 
store:  r:=p-i;  cdir:=o:  end 

ELSE  IF  (F7=0)  AND  (F6=255)  THEN  BEGIN  C:=C-1 
CDIR:=3:  ENO 

else  if  (f6=0)  and  (f5=255)  then  begin  r:=r*1 

store:  c:=c-i;  cdir:=2:  eno 
else  if  (f5=0)  and  (f4=255)  then  begin  r:=r+1 

CDIR:=2:  eno 
ELSE  IF  CF<.  =  0)  AND  (F3  =  255)  THEN  BEGIN  C:=C*1 

store;  r:=r*i;  cdir:=2:  eno 
ELSE  BEGIN  c:=c*i;  CDlR:=i:  end: 
enj: 
end: 
eno; 


PROCEDURE  INITIAL: 

BEGIN 

FOR  I  :=0  TO  255  DO  BEGIN 
CQLCID  :=o: 
ROUCI3  :=0! 

end: 

FOR  j:*0  TO  63  DO  BEGIN 
FOR  i:=0  TO  255  DO 
ACRtCD:=0: 

end; 
l:  =  o; 
end; 

procedure  first; 

VAR 

n,m  :integ;r; 

BEGIN 
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10-Dec-19&<.  13:53:33     VAX-11 
Source  Listing  6-Dec-198<>  17:38:13     _D.RA0: 

N  :=  o:  m:=o: 
FOR  C  :=0  TO  200  DO  BEGIN 
FOR  R  :=0  TO  63  DO  BEGIN 

IF  (FCR.C3=255)  UNO  CN  =  0)  THEN  BEGIN 

Ci  :=  C;  ri  :  =  r:  n  :=  n-m: 

END; 

IF  (FCR,255-C3=255)  AND  (M  =  0)  THEN  BEGIN 

C2  :  =  255-C:  R2  :=  R:  n    :=m*i; 
end: 
eno: 
end; 

wRITELNC"C0Ll=*,Cl,*R0Ul=',Ri,*C0L2«',C2,*R0*2=*.R2>: 

end: 
(#  main  program  *) 

9EGIN 

WRITELNC  'INPUT  CUT  OR  CONLINE  FILENAME*): 

REAOLN(NAME): 

OPEN  (INFILE, NAME, HISTORY  :=  OLD, 

ACCESS.METHOD  :  =  SE 3UENT I AL , 

RECORD_LENGTH  :  =  2 56  ,  RECORD_T YPE  :=FIXE0): 
OPEN  (OUTFILE.'TR. DAT', HISTORY  :  =  NE U , R ECORD_LENGTH  :=256, 

RECORD_TYPE  :=FIXED): 
RESETCINFILE) : 
REWRITE  (OUTFILEJ; 
R  :  =  0: 

WHILE  NOT  EOF  (INFILE)  00 
BEGIN 

READ  CINFILE.IMAGECR]); 

FOR  C  :=  0  TO  255  00 

fcr-»1.c*13  :=  imagecr,c3:        • 

r  :=  r*i  : 
end: 
first: 

R  :=Rl:  C  :  =  Ci:  CDIR  :  =  i: 
initial: 
while  (coc2)  do  begin 

store: 

cmove: 
end: 

FOR  l:=0  TO  COUNT  00 

acrowci3,c3lci33  :=255: 
for  r:=0  to  63  do  begin 
for  c:=c  to  255  do 

0utfilet.c3  :=  acr.cj: 
put<outfile): 
end: 

writelnc  * numb er=", count) ; 
close(in«=ile  ): 

ENO. 
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APPENDIX  C 
THE  SUBROUTINE  "PARAM" 
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SUBROUTINE  PARAMCX,T,Z,W,K,Z8,ZE,K,S.N,T,CX,CY,FP,IOPT,IPuR,IER) 
GIVEN  THE  SET  OF  OATA  POINTS  CXCD.YCI))  WITH  CORRESPONDING  Z- 

VALUES  ZCI)tI=lt2t H  AND  GIVEN  ALSO  THE  SET  OF  POSITIVE 

NUMBERS  WCI ) , 1=1 .2, ...M,  SUBROUTINE  PARAH  FINDS  A  SMOOTH  APPROXIMAT- 
ING CURVE  WITH  PARAMETER  REPRESENTATION  X  =  SX(Z),  Y  =  SYCZ). 
SXCZ)  AND  SYCZ)  ARE  TWO  SPLINE  FUNCTIONS  OF  DEGREE  K  WITH  THE  NUMBER 
ANC  THE  POSITION  OF  THE  KNOTS  T C J ) . J= 1 t 2 , . . . N  AUTOMATICALLY 
CHOSEN  3Y  THE  ROUTINE.  THE  SMOOTHNESS  OF  SX(Z)  AND  SYCZ)  IS 
ACHIEVED  BY  MINIMALIZING  THE  SUM(  DX  (  R  )*=*2  *0  Y  (  R  )S*2  )  WHERE  CXCR) 
AND  DY(R)  STANO  FOR  THE  D ISCONT I NUI T Y JUMP  OF  THE  KTH  DERIVATIVE 

OF  SXCZ)  ANO  SYCZ)  AT  THE  KNOT  TCR),  R=K*2. N-K-l. 

THE  AMOUNT  OF  SMOOTHNESS  IS  DETERMINED  BY  THE  CONOITION  THAT  FCP)  * 
SUMCWCI >*(( XCI)-SXCZCI)))**2  ♦CYCI)-SYCZCI)))**2))  BE  <  =  S.  WITH 
S  A  GIVEN  NON-NEGATIVE  CONSTANT. 
THE  SPLINE  FUNCTIONS  SXCZ)  AND  SYCZ)  ARE  GIVEN  IN  THEIR  B-SPLINE 

REPRESENTATION  C8-SPLINE  COEFFICIENTS  CXCJ),RESP.  CYCJ),J=1, N-K-l) 

AND  CAN  BE  EVALUATED  BY  MEANS  OF  FUNCTION  DERIV. 
CALLING  SEQUENCE: 

CALL  PARAMCX,Y,Z,W,M,Z9,ZE,K,S,N,T,CX,CY,FP,IGPT,IPAR,IER) 

INPUT  PARAMETERS: 

X      :  ARRAY, LENGTH  M,  CONTAINING  THE  ABSCISSAE  OF  THE  OATA  POINTS 
Y      :  ARRAY, LENGTH  M,  CONTAINING  THE  ORDINATES  OF  THE  OATA  POINTS 
W      :  ARRAY, HINIMUM  LENGTH  M, CONTAINING  THE  WEIGHTS  WCI). 
M      :  INTEGER  VALUE, CONTAINING  THE  NUMBER  OF  DATA  POINTS. 
K      :  INTEGER  VALUE, CONTAINING  THE  DEGREE  OF  SXCZ)  AND  SYCZ). 
S      :  REAL  VALUE, CONTAINING  THE  SMOOTHING  FACTOR. 
IOPT  :  INTEGER  VALUE  WHICH  TAKES  THE  VALUE  0  OR  1. 
IOPT=C:  THE  ROUTINE  WILL  RESTART  ALL  COMPUTATIONS. 
IOPT=l:  THE  ROUTINE  WILL  START  WITH  THE  KNOTS  FOUND  AT  THE 
LAST  CALL  OF  THE  ROUTINE.  IF  IOPT=l  THE  OUTPUT 
PARAMETERS  T  AND  N  ARE  INPUT  PARAMETERS  AS  WELL. 
IF  IOPT=l  THE  USER  MUST  PROVIDE  WITH  A  COMMON  BLOCK 
COMMON/ OPT1/NRDATACNE ST ),FPO,FPOLO,NPLUS 
INTEGER  FLAG. 


IPAR  : 
IPAR 


IPAR 

Z     : 

ZB.ZE: 


=  0:  FOR  EACH  DATA  POINT  CXCD.YCD)  THE  PROGRAM  AUTOMATICALLY 
CHOOSES  A  CORRESPONDING  VALUE  OF  THE  PARAMETER  Z,  I.E. 
ZC1)  =  0;ZCI)  =  ZCI-1>-»SQRTCCXCI)-XCI-1))*--*2«CYCI)-YCI-1))**2) 
THE  BOUNDARIES  FOR  THE  PARAMETER  Z  ARE  CHOSEN  AS  FOLLOWS 
ZB  =  ZC1)  :  ZE  =  ZCM). 

=  l:  THE  USER  HIMSELF  PROVIDES  WITH  THE  VALUES  OF  THE 
PARAMETER  Z  AND  WITH  THE  BOUNDARIES  ZB  AND  ZE. 

ARRAY, LENGTH  M,  CONTAINING  THE  VALUES  OF  THE  PARAMETER  Z 

CIPAR  =  1) 

REAL  VALUES, 

CIPAR  =  1). 


CONTAINING  THE  BOUNDARIES  OF  THE  PARAMETER  Z 


OUTPUT  PARAMETERS: 

T      :  ARRAY, LENGTH  NEST  CSEE  DATA 
WHICH  CONTAINS  THE  POSITION 
OF  THE  INTERIOR  KNOTS 
POSITION  OF  THE  KNOTS 


INITIALIZATION  STATEMENT), 

OF  THE  KNOTS, I.E.  THE  POSITION 

TCK+2), TCN-K-1),  AS  WELL  AS  THE 

TC  1)  =  TC2)  =  ...  =  TCK*1)  =  ZB  AND  ZE  = 


TCN-K)=. ..=TCN)  WHICH  ARE  NEEDEO  FOR  THE  B-SPLlNE  REPRESENT. 
CX.CY:  ARRAYS, LENGTH  NEST,  CONTAINING  THE  B-SPLlNE  COEFFICIENTS 

OF  SXCZ),  RESP.  SYCZ). 
N      :  INTEGER  V ALUE , CONT A I NI NG  THE  TOTAL  NUM3ER  OF  KNOTS. 
FP     :  REAL  VALUE,  WHICH  CONTAINS  THE  S UM C W I  * C X  I - S X C Z I ) >**2 ) 

♦  SUMCWKC  YI-SYCZI))**2),  1  =  1,2, ...M. 
IER    :  ERROR  COOE 

IER=0:  NORMAL  RETURN. 

I=R=-1 :NORMAL  RETURN, SXCZ)  ANO  SYCZ)  ARE  INTERPOLATING  SPLINES 
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IER=-2:N0RHAL  RETURN,  SX(2)  AND  STCZ)  ARE  POLYNOMIALS  OF  DEGREE  I". 
IER>0  :ABNORMAL  TERMINATION. 

I=R=l:THE  REQUIRED  STORAGE  SPACE  EXCEEDS  THE  AVAILABLE 

STORAGE  SPACE, SPECIFIED  BY  THE  PARAMETER  NEST. 

PROBABLY  CAUSEStNEST  OR  S  TOO  SMALL. 
IER=2:A  THEORETICALLY  IMPOSSIBLE  RESULT  WAS  FOUND  DURING 

THE  ITERATION  PROCESS. 

PROBABLY  CAUSESrTOL  TOO  SMALL 
IER=3:THE  MAXIMAL  NUMBER  OF  ITERATIONS  MAXIT  HAS  BEEN 

REACHEO. 

PROBABLY  CAUSESIMAXIT  OR  TOL  TOO  SMALL. 
IER=10:SOME  OF  THE  INPUT  OATA  ARE  INVALIOCSEE  RESTRICTIONS). 


.  M-l.    CIPAR  =  1) 


RESTRICTIONS: 

1)  M  >  K  >  0 

2)  IB    <=  Z(R)  <  Z(R-U)  <=  ZE,  R=l,2... 

3)  WCR)  >  0,  R=1,2,...M. 

4)  S  >=  0. 

5)  NEST  >=  2*K-»2. 

OTHER  SUBROUTINES  REQUIRED: 

eSPLlN.COSSIN, ROTATE, BACK, NKNOT, DISCO  AND  RATION. 
DIMENSION  X(M),Y(M),W(M),Z(M.),T(200},C.X(2CO).CYC200), 

<  FPINTC200).RX(200),RYC200),DIAGC200),DPRIMEC200), 

<  GC20C,6),8(200,7),QC40Q,6).H(7),NRDATA(200),A(200,5) 
C   COMMON/OPT 1/NRD AT  AC  NEST ),FPO ,FPOLO,NPLUS 

C      NRDATA:  INTEGER  ARRAY, LENGTH  NEST, WHICH  GIVES  THE  NUMBER  DF 

C  DATA  POINTS  INSIDE  EACH  KNOT  INTERVAL. 

C     FPO    :  REAL  VALUE,  WHICH  CONTAINS  THE  SUM(UI*( XI-SX (ZI ) 5**2) ♦ 

C  SUM(WI»(YI-SY(ZI))*S2)  WITH  SX(Z)  AND  SY(Z)  L E A S T- S 3U AR E S 

C  POLYNOMIALS  OF  DEGREE  K. 

C      FPOLD  :  REAL  VALUE, WHICH  CONTAINS  THE  SUM( W I* C XI- S X ( Z I >  )** 2  )♦ 

C  SUM(WI*(YI-SYCZI))**2)  WITH  SXCZ)  AND  SYCZ)  L E A S T- S3UAR E S 

C  SPLINE  FUNCTIONS  CORRESPONDING  TO  THE  LAST  FOUND  SET  OF 

C  KNOTS  BUT  ONE. 

C      NPLUS  :  INTEGER  VALUE, GIVING  THE  NUMBER  OF  KNOTS  OF  THE  LAST 

C  SET  MINUS  THE  NUMBER  OF  THE  LAST  SET  BUT  ONE. 

C0MM0N/0PT1/NRDATA, FPO, FPOLD, NPLUS 
C   DATA  INITIALIZATION  STATEMENT  TO  SPECIFY 

C      TOL    :  THE  REQUESTED  RELATIVE  ACCURACY  FOR  THE  ROOT  Oc  F(P)  =  S. 
C      MAXIT:  THE  MAXIMAL  NUMBER  OF  ITERATIONS  ALLOWED. 

C      NEST  :  AN  OVER-ESTIMATE  OF  THE  NUMBER  OF  KNOTS  N.  THIS  PARAMETER 
C  MUST  BE  SET  BY  THE  USER  TO  INDICATE  THE  STORAGE  SPACE 

C  AVAILABLE  TO  THE  SUBROUTINE.  THE  DIMENSION  SPECIFICATIONS 

C  OF  THE  ARRAYS   T , C X , C Y , NRD A T A , F P I N T , R X , R Y , D I  AG , 0  PR  I  ME ( N) , 

C  A(N,K) ,G(N,K*1 ),B(N,K*2) ,SCM,K*1 )  AND  H(K*2)  DEPENO 

C  ON  N,M  AND  K.  SINCE  N  IS  UNKNOWN  AT  THE  TIME  THE 

C  USER  SETS  UP  THE  DIMENSION  INFORMATION  AN  OVER-ESTIMATE 

C  OF  THESE  ARRAYS  WILL  GENERALLY  BE  MADE.  THE  FOLLOWING 

C  REMARKS  ARE  INTENDED  TO  HELP  THE  USER 

C  1)  2*K-»2  <=  N  <=  M*K*I 

C  2)  THE  SMALLER  THE  VALUE  OF  S,  THE  GREATER  N  WILL  BE. 

C  3)  NORMALLY  N  =  M/2  IS  AN  OVER-ESTIMATE. 

DATA  TOL/0.001/iMAXIT/20/, NEST/200/ 
C    BEFORE  STARTING  COMPUTATIONS  A  DATA  CHECK  IS  MADE.  IF  THE  INPUT 
C   DATA  ARE  INVALID  CONTROLE  IS  IMMEDIATELY  REPASSED  TO  THE  ORIVER 
C   PROGRAM  (IER=10). 

IER  =  0 

KI  =  K+l 

K2  *  Kl*l 

NMIN  =  2*K1 


1540 
1550 
1560 
1570 
1530 
1590 
1600 
1610 
1620 
1630 
1640 
1650 
1660 
1670 
1680 
1690 
1700 

mo 

1720 
1730 
1740 
1750 


1790 
1600 
1810 
1820 
1830 
1840 
1650 
1860 
1870 
1830 
1P90 
1900 
1910 
1920 
1930 
1°40 
1950 
1960 
1970 
1980 
1990 
2000 
2010 
2020 
2030 
2040 
2050 
2060 

2080 
2090 
2100 
2110 
2120 
2130 
?1  4C 


119 


ERPOLATING  SPLINES. 

N  =  NMAX 
ST  WHETHER  THE  REQUIRED  STORAGE  SPACE  EXCEEDS  THE  AVAILABLE  ONE. 

IF(N.GT.NEST)  GO  TO  <*20 
NO  THE  POSITION  OF  THE  INTERIOR  KNOTS  IN  CASE  OF  INTERPOLATION. 

HM  =  M-Kl 


TEST 
FI 


0  60 


30 


20 
30 


MM  =  M-M 
IFCMK1.EQ.0)  GO  T 
K3  =  K/2 
I  =  K2 
J  =  K3+2 

IFCK3*2.EC.K )  GO  TO 
DO  20  L=1,MK1 

TCI)  =  ZCJ  ) 

I  =  1*1 

J  =  J*  1 
CONTINUE 
GO  TO  60 
DO  40  L=1,MK1 

TO  =  (Z(J)-Z(J-1))*0, 

I  =  1*1 

J  =  J*l 
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FCZ) 


80 


C   FETCH 


DIAGCI)  =  0 

RXCI)  =  0. 

RYCI  )  =  0. 

DO  80  J  =  1,K. 
ACI.J)  =  0. 
CONTINUE 
L  =  M 
00  130 


0  130  IT=1,M 

THE  CURRENT  DATA  POINT  X C I T  )  ,  Y (  I  T  )  ,  Z C  I  T  )  . 

XI  =  XCIT) 

YI  =  YCIT) 
Z(IT) 


C   EVA 


ZI  =  Z(IT) 
WI  =  W(IT) 
iRCH  FOR  KNOT  INTERVAL  T(L)  <  =  ZI  <=  T(L*1). 
IFCZI.GE.TCL  +  1)  .AND.  L.NE.NK.1)  L  =  L*l 
iLUATE  THE  CK-H)  NON-ZERO  8-SPLINES  AT  ZI  AND  STORE  T 
CALL  BSPLIN(T,N,H,ZI.L,H) 
DO  90  1=1,  Kl 

IF(HCI).LT.0.1E-07)  H(I)  =  0. 
OC  IT,  I  )  =  HCI) 
0        CONTINUE 
ROTATE  THE  NEW  ROW  OF 
GIVENS  TRANSFORMATION 
J  =  L-Kl 


HEM  IN  Q. 


J 

DO 


L-M 
110  1=1 ,K1 


THE  OBSERVATION  MATRIX  INTO  TRIANGLE 
S  WITHOUT  SQUARE  ROOTS. 


BY 


110  1=1 ,K1 
IFCWI.EO.O.)  GO  TO  130 
J  =  J*l 


PIV  =  H(I) 
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IFCPIV.EQ-O.)  GO  TO  110 
C   CALCULATE  THE  PARAMETERS  OF  THE  GIVENS  TRANSFORMATION. 

CALL  COSSIN(PIV,WI,DIAG(J),COS,SIN) 
C   TRANSFORMATIONS  TO  RISHT  HAND  SIDES. 

CALL  ROTATE(PIV,COS,SIN,XI,RX(J)> 

CALL  ROTATECPIV,COS,SIN,YI,RY(J)) 

IFCI.EQ.K1)  GO  TO  120 

12  =  0 

13  =  1*1 

DO  100  II  =  13, Kl 
12  =  12*1 
C   TRANSFORMATIONS  TO  LEFT  HAND  SIDE. 

CALL  ROTATECPIV,COS.SIN,H<Il),A(J,I2)) 
100  CONTINUE 

110        CONTINUE 
C   ADD  CONTRIBUTION  OF  THIS  ROW  TO  THE  SUM  OF  SQUARES  OF  RESIOJAL 
C   RIGHT  HAND  SIDES. 
120        FP  =  FP+WI*<XI**2*TI**2> 
130     CONTINUE 

IFCIER.EQ.-2)  FPO  =  FP 
C   BACKWARD  SU3STITUTION  TO  OBTAIN  THE  B-SPLINE  COEFFICIENTS. 
CALL  3ACK(A,RX,NK1,K,CX) 
CALL  3ACK(A,RY,NK1,K,CY) 
C   TEST  WHETHER  THE  APPROXIMATION  X = SXI NF ( Z ) , Y=S Y INF C Z )  IS  AN 
C   ACCEPTABLE  SOLUTION. 
FPMS  =  FP-S 

IF(ABS(FPMS).LT.ACC)  GO  TO  44 0 
C   IF  F(P=INF)  <  S  ACCEPT  THE  CHOICE  OF  KNOTS. 

IFCFPMS.LT.O.)  GO  TO  250 
C   IF  N=NMAX.SXINF(Z)  AND  SYINF(Z)  ARE  INTERPOLATING  SPLINES. 

IF(N.EQ-NMAX)  GO  TO  430 
C   INCREASE  THE  NUMBER  OF  KNOTS. 

C   IF  N=NEST  WE  CANNOT  INCREASE  THE  NUMBER  OF  KNOTS  BECAUSE  OF 
C   THE  STORAGE  CAPACITY  LIMITATION. 

IF(N.EQ.NEST)  GO  TO  420 
C   DETERMINE  THE  NUMBER  OF  KNOTS  NPLUS  WE  ARE  GOING  TO  AOD. 
IFCIER.E3.0)  GO  TO  140 
NPLUS  =  1 
IER  =  0 
GO  TO  150 
140      NPL1  =  NPLJS*2 

IF(FPOLD-FP.GT.ACC)  NPL1  =  FLO AT( NPLUS )*FPMS/< F POLD-F P ) 
NPLUS  =  MIN0CNPLUS*2,MAX0(NPL1,NPLUS/2,D) 
150     FPOLD  =  FP 
C   COMPUTE  THE  SUM(WI*(CXI-SXINF(ZI))**2*CYI-SYINFCZI)>**2))  FOR 
C   EACH  KNOT  INTERVAL  T(J*K)  <=  ZI  <=  T  (  J  ♦  K  ♦  1  )  ANO  STORE  IT  IN 

C   FPINTCJ),J=1,2, NRINT. 

FPART  =  0. 

1  =  1 

L  =  K2 

NEW  =  0 

00  180  IT=1,M 

IF(ZCIT).LT.T(L;  .OR.  L.GT.NKD  GO  TO  160 
NEW  =  1 
L  =  L*l 
160        TERM1  =  0. 
TERM2  =  0. 
LO  =  L-K2 
DO  170  J=1,K1 
LO  =  L0*1 
TERMl  =  TERMl*CX(LO)sCCIT, J) 
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TERM2  =  TERM2*CYCL0)*CCIT, J) 
170        CONTINUE 

TERM  =  W(IT)*((TERMl-xCIT))$?2*( TERM2-YC IT ))s*2) 
FPART  =  FPART+TERM 
IFCNEW.E3.0)  GO  TO  180 
STORE  =  TERM*0.5 
FPINTCI)  =  FPART-STORE 
I  =  1*1 
FPART  =  STORE 
NEW  =  0 
180      CONTINUE 

FPINT(NRINT)  =  FPART 
00  190  L=1,NPLUS 
C   AOO  A  NEw  KNOT. 

CALL  NKN3T(Z,M,T,N,FPINT,NR0ATA,NRINT) 
C   TEST  WHETHER  HE  CANNOT  FURTHER  INCREASE  THE  NUM3ER  OF  KNOTS. 
IFCN.E3.NMAX  .OR.  N.EQ.NEST)  GO  TO  200 
190      CONTINUE 
C   RESTART  THE  COMPUTATIONS  WITH  THE  NEW  SET  OF  KNOTS. 

200   CONTINUE 
C   TEST  WHETHER  THE  APPROXIMATION  X  =  S X C Z  )  t T= SY ( Z )  WITH  SX(Z)  AND  ST(Z) 
C   THE  LEAST-SQUARES  KTH  DEGREE  POLYNOMIALS,  IS  A  SOLUTION  OF  OUR  PROBLEM 

250   IFCIER.E3.-2)  GO  TO  440 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C   PART  2:  DETERMINATION  OF  THE  SMOOTHING  SPLINES  SXP(Z)  AND  SYP(Z).    C 

C   WE  HAVE  DETERMINED  THE  NUMBER  OF  KNOTS  AND  THEIR  POSITION.  C 

C  WE  NOW  COMPUTE  THE  B-SPLINE  COEFFICIENTS  OF  THE  SMOOTHING  SPLINES  C 
C  SXP(Z)  ANO  SYP(Z).  THE  OBSERVATION  MATRIX  A  IS  EXTENDED  BY  THE  ROWS  C 
C   OF  MATRIX  B  EXPRESSING  THAT  THE  KTH  DERIVATIVE  DISCONTINUITIES  OF     C 

C   SXP(Z)  AND  SYP(Z)  AT  THE  INTERIOR  KNOTS  T(K*2), TCN-K-1)  MUST  BE   C 

C  ZERO.  THE  CORRESPONDING  WEIGHTS  OF  THESE  AOOITIONAL  ROWS  ARE  SET  C 
C  TO  1/SORTCP).  ITERAT1VELY  WE  THEN  HAVE  TO  DETERMINE  THE  VALUE  OF  P  C 
C  SUCH  THAT  F(P)=SUM(WI*((XI-SXP(ZI))S*2*CYI-SYP(ZI))**2)  BE  =  S.  WE  C 
C  ALREADY  KNOW  THAT  THE  L c A  ST- SOU  ARE S  POLYNOMIALS  CORRESPOND  TO  P  =  0,  C 
C  AND  THAT  THE  L E A  ST- S QU A R E S  SPLINES  CORRESPOND  TO  P=INFINITY.  THE  C 
C  ITERATION  PROCESS  WHICH  IS  PROPOSED  HERE,  MAKES  USE  OF  RATIONAL  C 
C    INTERPOLATION.  SlNCS  FCP)  IS  A  CONVEX  ANO  STRICTLY  DECREASING  C 

C   FUNCTION  OF  P,  IT  CAN  BE  APPROXIMATED  BY  A  RATIONAL  FUNCTION  C 

C  R(P)  =  (U*P»V>/(P*W).  THREE  VALUES  OF  P<P1,P2,P3)  WITH  CORRESPOND-  C 
C  ING  VALUES  OF  F(P)  CF1=F(P1)-S,F2=F(P2)-S,F3=P(&3)-S)  ARE  USED  C 
C  TO  CALCULATE  THE  NEW  VALUE  OF  P  SUCH  THAT  R(P)=S.  CONVERGENCE  IS  C 
C   GUARANTEED  BY  TAKING  F1>0  AND  F3<0.  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C    EVALUATE  THE  DISCONTINUITY  JUMP  OF  THE  KTH  DERIVATIVE  OF  THE 
C   B-SPLINES  AT  THE  KNOTS  T(L),L=K*2, N-K-l  ANO  STORE  IN  B. 

CALL  DISCOCT.N, K2,8) 
C   INITIAL  VALUE  FOR  P. 

PI  =  0. 

Fl  =  FPO-S 

P3  =  -1. 

F3  *  FPMS 

P  =  -F1/F3 

ICHECK  =  0 

N8  =  N-NMIN 
C   ITERATION  PROCESS  TO  FIND  THE  ROOT  OF  F(P)  =  S. 

DO  350  ITER=1,MAXIT 
C   THE  ROWS  OF  MATRIX  S  WITH  WEIGHT  1/SORTCP)  ARE  ROTATED  INTO 
C   THE  TRIANGULARISED  OBSERVATION  MATRIX  A  WHICH  IS  STORED  IN  G. 
PINV  =  1.0/P 
DO  260  1=1, NK1 
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DPRIMECI)  =  0I4GCI)  4590 

CXCI )  =  RXCI  )  4600 

CYCI)  =  RYU)  4610 

GCI.M)  =  0.  4620 

DO  260  J=1,K  4630 

G(I ,J)  =  A(I,J)  4640 

260     CONTINUE  4650 

DO  300  IT=1,N8  4660 

C   THE  ROW  OP  MATRIX  B  IS  ROTATED  INTO  TRIANGLE  BY  GIVENS  TRANSFORMATIONS.   4670 

DO  270  1=1, K2  4680 

H(I)  =  BCIT.I)  4690 

270        CONTINUE  4700 

XI  =  0.  4710 

YI  =  0.  4720 

Ul  =  PINV  4730 

DO  290  J=IT,NK1  4740 

IFCWI.EQ.O.)  GO  TO  300  4750 

PIV  =  H(l)  4760 

C   CALCULATE  THE  PARAMETERS  OF  THE  GIVENS  TRANSFORMATION.  4770 

CALL  COSSINCPIV.UI.DPRIMEC J), COS, SIN)  4780 

C   TRANSFORMATIONS  TO  RIGHT  HAND  SIDES.  4790 

CALL  ROTATECPIV,COS,SIN,XI,CX( J))  4800 

CALL  R3TATECPIV,C0S,SIN,YI,CYCJ)>  4810 

IF(J.EO.NKl)  GO  TO  300  4820 

12  =  K.1  4330 

IFCJ.GT.N8)  12  =  NK1-J  4840 

DO  280  1=1,12  4850 

C   TRANSFORMATIONS  TO  LEFT  HANO  SIDE.  4860 

CALL  R3TATE(PlV,C0S,SlN,H(I-fl),G(J,I))  4870 

rt(I)  =  H(I*1)  4830 

280          CONTINUE  4890 

HCI2+1)  =  0.  4900 

290        CONTINUE  4910 

300      CONTINUE  4920 

C   BACKWARO  SUBSTITUTION  TO  OBTAIN  THE  B-SPLINE  COEFFICIENTS.  4930 

CALL  BACH(G,CX,NK1,K,CX)  4940 

CALL  BACK(G,Ct,NK,l,K,CY)  4950 

C    COMPUTATION  OF  F(P).  4960 

FP  =  0.  4970 

L  =  K2  4980 

00  330  IT=1,M  4990 

IFCZCIT).LT.T(L)  .OR.  L.GT.NK1)  GO  TO  310  5030 

L  =  L-H  5010 

310        LO  =  L-K 2  5  02  0 

TERM1  =  0.  5030 

TERM2  =  0.  5040 

DO  320  J=1,K1  5050 

LO  =  L0*1  5060 

TERM1  =  TERM1*CX(LO)*OCIT, J)  5070 

TERM2  =  TERM2-»CY(L0)«Q(IT,J)  5090 

320        CONTINUE  5090 

FP  =  FP*-W(IT)*(CTERM1-X(IT))«2-»(TERM2-YCIT)  )**2)  5100 

330      CONTINUE  5110 

C   TEST  WHETHER  THE  APPROXIMATION  X  =  SXP C Z ) , Y  =  S YP < 1  )  IS  AN  ACCEPTABLE  5120 

C   SOLUTION.  5130 

FPMS  =  FP-S  5140 

IFCABS(FPMS).LT.ACC)  GO  TO  440  5150 

C   TEST  WHETHER  THE  MAXIMAL  NUMBER  OF  ITERATIONS  IS  REACHED.  5160 

IF(ITER.EQ.MAXIT)  GO  T3  400  ■  5170 

C   CARRY  OUT  ONE  MORE  STEP  Or     THE  ITERATION  PROCESS.  5130 

P2  =  P  5  190 
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F2  =  FPMS 

IFC  ICHEC/C.NE.O)  GO  TO  340 

IFCCF2-F3).GT.ACC)  GO  TO  335 
OUR  INITIAL  CHOICE  OF  P 

P  =  P*0.1E-02 

P3  =  P2 

*3  =  F2 

GO  TO  350 
335     IFCCF1-F2).GT.ACC)  GO  TO  3*0 
OUR  INITIAL  CHOICE  OF 

TYPE  *,  'VALUE  OF 

P  =  P*0.1E*04 

PI  =  PZ 


15  TOO  LASGE. 


P  IS  TOO  SHALL 
P-.P 


Fl  =  F2 
GO  TO  350 
C   TEST  WHETHER  THE  ITERATION  PROCESS  PROCEEDS  AS  THEORETICALLY 
C   EXPECTED. 
340     IFCF2.GE.F1  .OR.  F2.L5.F3)  GO  TO  410 
ICHECK  =  1 
C   FIND  THE  NEW  VALUE  FOR  P. 

P  =  RATIONCPl.Fl, P2.F2.P3.F3) 
350   CONTINUE 
C   ERROR  COOES  AND  MESSAGES. 
400   IER  =  3 

GO  TO  440 
2 


410 

420 

430 
440 


IER  = 
GO  TO  440 
IER  =  1 
GO  TO  440 
IER  =  -1 
RETURN 
ENO 

SUBROUTINE  B SPLI NC T , N, K , X , L , H ) 
C   SUBROUTINE  BSPLIN  EVALUATES  THE  C**l>  NON-ZERO  B-SPLINES  OF 
C   DEGREE  K  AT  TCL)  <=  X  <  TCL*1)  USING  THE  STABLE  RECURRENCE 
C   RELATION  OF  OE  BOOR  ANO  COX. 

C   THE  DIMENSION  SPECIFICATIONS  OF  THE  FOLLOWING  ARRAYS  MUST  BE 
C   AT  LEAST  HC  K  ♦  1  )  .  MM  C  C.  )  . 

DIMENSION  TCN),HC6),HHC5) 
HC1)  =  1. 
00  20  J=1,K 
DO  10  1  =  1, J 

HHCI)  =  HCI) 
10     CONTINUE 
HC1)  =  0. 
DO  20  1  =  1, J 
LI  =  L*I 
LJ  =  LI-J 

F  =  HHCD/CTCLI  )-TCLJ)) 
H(I)  =  HCI)*F»CTCLI)-X) 
MCI*1)  =  FSCX-TCLJ)) 
20   CONTINUE 
RETURN 
ENO 

SUBROUTINE  COSSINCPIV.WI.WW.COS.SIN) 
C   SUBROUTINE  COSSIN  CALCULATES  THE  PARAMETERS  OF  A  GIVENS 
C   TRANSFORMATION  WITHOUT  SQUARE  ROOTS. 
STORE  =  PIVOWI 
DO  =  WW*ST0RE*PIV 
IFCAS5C0D3.LT. l.E-36)  DD=l.E-36 
CCS  =  WW/DO 
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SOLUTION  OF  THE  SYSTEM  OF 

N  UNIT  UPPER  TRIANGULAR  MATRIX 


SIN  =  STORE/OO 

uw  =  DO 

HI  =  COSSWI 

RETURN 

ENO 

SUBROUTINE  ROTATE<PIV,COS,SIN,A,3> 
SUBROUTINE  ROTATE  APPLIES  A  GIVENS  ROTATION  TO  A  ANC  B, 

STORE  =  B 

B  =  COS*STORE*SIN*A 

A  =  A-PIV*STOR£ 

RETURN 

END 

SUBROUTINE  3ACKCA,2,N,K,C> 
SUBROUTINE  BACK.  CALCULATES  THE 
EQUATIONS  A$C  =  2  WITH  A  A  N  X 
OF  BANDWIDTH  K-»l. 

ATTENTION:  THE  FIRST  OIMENSION  SPECIFICATION  OF  MATRIX  A  MUST 
BE  THE  SAME  AS  IN  THE  CALLING  PROGRAM. 

DIMENSION  A(200,K),2(N),C(N) 

C(N)  =  2CN) 

I  =  N-l 

IFCI.EQ.O)  GO  TO  30 

00  20  J=2,N 
STORE  =  2(1) 
II  =  K 

IFCJ.LE.K)  II  *  J-l 
M  =  I 

00  10  L=l ,11 
M  =  M*l 

STORE  =  STORE-C(M)SA(I,L) 
CONTINUE 
CCI)  =  STORE 

1  =  1-1 
CONTINUE 
RETURN 
ENO 

SUBROUTINE  NKNOT(X,M,T,N,FPINT,NROATA,NRINT) 
SUBROUTINE  NKNOT  LOCATES  AN  ADDITIONAL  KNOT  FOR  A 
K  AND  ADJUSTS  THE  CORRESPONDING  PAR AMET ER S, I . c. 
T       :  THE  POSITION  OF  THE  KNOTS. 
N       :  THE  NUMBER  OF  KNOTS. 
NR1NT  :  THE  NUMBER  OF  K NO T I  NT E R V  A L S . 
FPINT  :  THE  SUM  OF  SQUARES  OF  RESIDUAL  RIGHT  HAND 

FOR  EACH  KNOT  INTERVAL. 
NRDATA:  THE  NUMBER  OF  DATA  POINTS  INSIDE  EACH  KNOT  INTERVAL. 
THE  ARRAYS  T.FPINT  ANO  NRDATA  MUST  HAVE  THE  SAME  DIMENSION 
SPECIFICATIONS  AS  IN  THE  CALLING  PROGRAM. 

OIMENSION  X(N},T(200),FPINT(200),NRDATA(200> 
K  =  (N-NRINT-1 )/2 
SEARCH  FOR  KNOT  INTERVAL  T(NUMBER*K)  <=  X 
FPINTCNUMBER)  IS  MAXIMAL  ON  THE  CONOITION 
NOT  EQUALS  2ER0. 
FPMAX  =  0. 
JBEGIN  =  1 
DO  20  J=1,NRINT 

JP01NT  =  NRDATA(J) 
IFCFPMAX.GE.FPINTCJ) 
FPMAX  =  FPINTCJ) 
NUMBER  =  J 
MAXPT  =  JPOINT 
MAXBEG  =  J3EGIN 


10 


20 
30 


SPLINE  OF  DEGRF.E 


SIDES 


<=  T(NUMBER*K*1 )  WHERE 
THAT  NRDATA(NUMBER) 


.OR.  JPOINT. EQ.O)  GO  TO  10 
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10      JBEGIN  =  JBEGIN*JPOINT*l 
20   CONTINUE 
C   LET  COINCIDE  THE  NEW  KNOT  T ( NUM3E R *K *1 )     WITH  a  DATA  PJINT  X(NRX) 
C   INSIDE  THE  OLD  KNOT  INTERVAL  T(NUMBER*X)  <=  X  <=  T ( NUMB E R ♦ K ♦  1 ) . 
IHALF  =  MAXPT/2*1 
NRX  =  MAXBEG+IHALF 
NEXT  =  NUMBER+1 
IF(NEXT.GT.NRINT)  GO  TO  40 
C   A0JU5T5  THE  DIFFERENT  PARAMETERS. 
00  30  J=NEXT,NRINT 
JJ  =  NEXT+NRINT-J 
FPINT(JJ»1)  =  FPINT(JJ) 
NRDATACJJ+1)  *  NROATA(JJ) 
JK  =  JJ»X 
T(JK*1)  =  TCJK) 
30   CONTINUE 
40   NRDATA(NUMBER)  =  IMALF-1 

NROATA(NEXT)  =  MAXPT-IHALF 

FPINTC NUMBER)  =  FPMA X*FL OAT ( NROA T A ( NUMB E R ) )/FLO AT ( M A XP T ) 
FPINT(NEXT)  =  FPMAX*FL0AT(NROATA(NEXT) )/FL 0 A T ( MA  X P T ) 
JK  =  NEXT*K 
TCJK)  =  X(NRX) 
N  =  N*l 

NRINT  =  NRINT*1 
RETURN 
END 

SUBROUTINE  DISC0(T,N,K2,B) 
C   SUBROUTINE  DISCO  CALCULATES  THE  DISCONTINUITY  JUMPS  Oc  THE  KTH 
C   DERIVATIVE  OF  THE  B-SPLINES  OF  DEGREE  K  AT  THE  KNOTS  T ( X ♦ 2 )  .  - T ( N-K- 
C   THE  FIRST  DIMENSION  S ° E C I FI C A T I  ON  OF  THF  MATRIX  E  MUST  3E  THE  SAME 
C    IN  THE  CALLING  PROGRAM;  H  MUST  HAVE  A  DIMENSION  SPECIFICATION  AT 
C   LEAST  2*K*2. 

DIMENSION  T(N),B(200,K2),M(12) 
XI  =  K2-1 
K  =  Kl-1 
NR1  =  N-Xl 
DO  40  L=K2,NK1 
LMK  =  L-Kl 
00  10  J=1,X1 
IX  =  J*K1 
LJ  =  L*J 
LX  =  LJ-K.2 
H(J)  =  TCL)-T(LK) 
H(IK)  =  T(L)-T(LJ) 
10     CONTINUE 
LP  =  LMK 
DO  30  J=1,K2 
JX  =  J  +  K 
PROO  =  1. 
DO  20  I=J,JK 

PROD  *  PROD*H(I) 
20       CONTINUE 

LK  r  LP«-K1 

B(LHK,J>  *  (T(LK)-T(LP))/PROD 
LP  =  LP*1 
30     CONTINUE 
40   CONTINUE 


1) 
AS 


RETURN 
END 

FUNCTION 
GIVEN  THREE 


RATION (PI  ,  F1,P2,F2,P3,F3) 

POINTS  (P1.F1), (P2.F2)  AND  (P3.F3),  FUNCTION  RATION 
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GIVES  THE  VALUE  OF  P  SUCH  THAT  THE  RATIONAL  INTERPOLATING  FUNCTION 
OF  THE  FORM  R(P)  =  (U3P*V)/(P-»-W)  EQUALS  ZERO  AT  P. 

IFCP3.GT.0.)  GO  TO  10 
VALUE  OF  P  IN  CASE  P3  =  INFINITY. 

P  =  CP1*(F1-F3)*F2-P2*CF2-F3>*F1)/C(F1-F2)*F3) 

GO  TO  20 

P3  "=  INFINITY. 


VALUE 

10   HI 

H2 

H3 


OF  P  IN  CASE 
=  F1$(F2-F3) 
=  F2*(F3-F1) 
=  F3«(F1-F2) 


P  =  -(  P13P2*H3-»P2<-P3*H1*P3=»P1*H2  )/CPl*Hl  *P2*H2»P3*H3  ) 


AOJUST  THE  VALUE  OF  P1.F1.P3  AND  F3  SUCH  THAT  Fl 
20   IFCF2.LT.0.)  GO  TO  30 


>  0  ANO  F3  <  0. 


30 


40 


PI 
Fl 
GO 
P3 
F3 


C 

C 
C 

c 
c 
c 
c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


=  P2 
=  F2 
TO  40 
=  P2 
=  F2 

RATION  = 

RETURN 

ENO 

FUNCTION 


OERIVCT.N.C.NKl.NU.ARG.L) 


DEGREE  K 
OR 


FUNCTION  OERIV  EVALUATES  A  SPLINE  S(X)  OF 
GIVEN  IN  ITS  NORMALIZED  B-SPLlNE  REPRESENTATION 
DERIVATIVES  OF  ANY  SPECIFIED  ORDER  NU. 

CALLING  SEQUENCE 

VALUE  =  OERI  V(T,N,C,NU,  NU.ARG.L) 


WHICH  IS 
CALCULATES 


INPUT 
T 


N 

C 

NK1 

NU 

ARG 

L 


CONTAINS 
POSITION 


THE  POSITION 
OP  THE  INTERIOR 
AS  THE  POSITION  OF  THE 
..TCN)  WHICH  ARE  NEEDED 


PARAME  TERS : 
:  ARRAY .MINIMUM  LENGTH  N,  WHICH 
OF  THE  KNOTS  OF  S(X),I.E.  THE 
KNOTS  T(K*2),...T(N-K-1)  AS  WELL 
KNOTS  T(l  )  ,...T(K-»1)  ANO  T(N-K},. 
FOR  THE  B-SPLINE  REPRESENTATION. 

INTEGER  VALUE  GIVING  THE  TOTAL  NUMBER  OF  KNOTS  OF  S(X). 
ARRAY, LENGTH  NK1,  CONTAINING  THE  B-S°LlNE  COEFFICIENTS. 
INTEGER  VALUE, GIVING  THE  DIMENSION  OF  S(X),I.E.NK1  =  N-K-l. 
INTEGER  VALUE  WHICH  GIVES  THE  ORDER  OF  THE  DERIVATIVE. 
REAL  VALUE. GIVING  THE  VALUE  OF  THE  ARGUMENT. 

INTEGER  VALUE  WHICH  SPECIFIES  THE  POSITION  OF  THE  ARGUMENT 
I.E.    TCL)  <=  ARG  <   T(L*1)    OR 
L  =  NK1  IF  ARG  =  T(NKl-fl). 


OUTPUT  PARAMETER: 

VALUE:  REAL  VALUE. GIVING 
SCX)  AT  X  =  ARG. 


THE  VALUE  OF  THE  NUTH  DERIVATIVE  OF 


OTHER  SUBROUTINES  REQUIRED:  NONE. 

RESTRICTIONS: 

1)  NU  >=  0 

2)  TCK«1)  <=  ARG  <=  T(NK1*1) 

THE  DIMENSION  SPECIFICATION  OF  THE  AR3AY  H  MUST  BE  AT  LEAST  K*l. 
OIMENSION  T(N),CCNK1),H(6) 
OERIV  =  0. 
Kl  *  N-NK1 

IFCNU.LT.O  .OR.  NU.GE.K1)  RETURN 
DO  100  1=1, Kl 
IK  =  L*I-K1 

H(I)  =   CCIK) 
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128 


100 


200 


300 


400 
500 


600 


=  CMCI)-HCI-1))/CTCLJ)-TCLD) 


IU-1)  GO  TO  500 


CONTINUE 

IFCNU.EQ.O)  GO  TO  300 
NU1  =  NU*1 
00  200  J=2,NU1 
00  200  JJ=J,K1 

I  -     J-»M-JJ 

Li  =  L+I-Kl 

LJ  =  L*I-J*1 

H(I) 
CONTINUE 
IFCNU-EG. 
NU2  =  NU  +  2 
DC  400  J=NU2,K1 
DO  400  JJ=J,K1 

I  =  J+K1-JJ 

LI  =  L*I-K1 

LJ  =  L*I-J*1 

H(I)  =  C(AR&-T(Ll))*HCI)*(T(LJ)-ftRG)*MCI-l))/CT(LJ)-TCLI)) 
CONTINUE 
DERI  V  =  HCM  ) 
IF(NU.EC.O)  RETURN 
DO  600  1=1, NU 

DERIV  =  OERIV*FLOAT(M-I> 
CONTINUE 
RETURN 
END 
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APPENDIX  D 
THE  PROGRAM  TO  FIND  B-SPLINE  COEFFICIENT 
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Source  Listing 

PROGRAM  3SPLINECINPUT, OUTPUT, I NFILE.OUTFILE): 


10-Dec-193<.  16:37:51 
10-Dec-1984  16:37:39 


VAX-1 
DRAO 


TYPE 

BYTE  =  0..255: 

IMAGEROWl  =  PACKED  ARRAY  CO. .2553  OF  BYTE: 

ROW1  =  PACKED  ARRAY  CO. .2573  OF  BYTE; 

SHIP  =  PACKED  ARRAY  C C . . 63 , 0 . . 25 5 3  OF  BYTE; 

DAI  =  PACKED  ARRAY  CO. .5123  OF  REAL: 

DA2  =  PACKED  ARRAY  CI. .3003  OF  REAl; 
VAR 

R,C,F0,F1,F2,F3,F<.,F5,F6,F7,F8   :   BYTE; 

F  :  ARRAY  CO. .653  OF  ROWl; 

A, IMAGE  :  SHIP; 

INFILE  :  FILE  OF  IMAGEROWl: 

SPX.SPY, SPX1.SPY1  :  PACKED  ARRAY  CO. .5123  OF  REAL! 

R1,C1,R2,C2,CDIR,I,J,J1,C0UNT,NEG  : INTEGER; 

TEMPX.TEHPY ,RA  :  REAL: 

M, I0PT.K ,IPAR,N,IER,ANS,NU,ANS1 ,ANS2, ANS3  :  INTEGER: 

NK1 .NENO.L.MET  :INTEGER; 

u,;  :packed  array  cc. 5123  of  real; 

S.ZB.ZE.FP  :  real; 

ARG, THETA.TOLE  :  REAL: 

FLAG_BEGIN,PLAG_ENO,AK,LUMP  :  INTEGER; 

BEG, In  :  PACKED  ARRAY  CI. .53  OF  REAL: 

BEGN.ENN  :  PACKED  URRAY  CI. .53  OF  INTEGER; 

T,CX,CY  :PACKED  ARRAY  CI. .3003  OF  REAL! 

COL, ROW  :PACKED  ARRAY  CO. .5123  OF  REAL: 

X,Y  :PACKED  ARRAY  CO.. 5123  OF  REAL: 

DCY.OCX.ARE ,CY1 .DMAX.DMIN  :  REAL: 

CYMIN.CYMAX ,MAxCY  :  PACKED  ARRAY  CI. .53  OF  REAL: 

AREA  :  PACKED  ARRAY  CI.. 53  OF  REAL: 

OUTFILE  :  FILE  OF  IMAGEROWl: 

NAME  :  PACKEO  ARRAY  CI. .203  OF  CHAR; 

PEAK, STAR, TER  :PACKED  ARRAY  CI. .53  OF  REAL: 

(*  FILTER  THE  POINTS  ft) 


PROCEOURE  STORE: 
BEGIN 

TEM 

IF 

F 


C 
R 
M 
I 
5NO 
ELS 
C 
R 

M 
I 

ENO 

end; 


Px: 

I>0 

OR 
IF 
TH 

no; 

OLC 
Owe 
:  =  I 
:  =  I 


=c-i:  tempy:  =  6<.-cr-i); 

THEN  BEGIN 
j:=0  TO  M  DO  EEGIN 
CCOLCJ3=TEMPX)  AND  ( ROW C J 3=TEMP Y ) 
EN  i:=j: 


13  :=TEmpx; 
13  :=tempy; 

♦  l : 


E  BEGIN 

DLCI3  :=TEMPX: 


OWCI3 
:  =  0: 

:  =  i: 


=tempy; 
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iourc  e  Listing 


10-Dec-198<.  16:37:51 
10-Dec-1984  16:37:39 


VAX-1 

_drao 


PROCEDURE  CMOVE; 
BEGIN 

FO  :=F[R-1.C3:  Fl:  =  FCR-l,C  +  l3:  F2  :  =  FCR . C ♦ 1 3 : 
F3:  =  FCR*1,C*13:  F4:=FCR*1  .C3;  F5:=FCR*1 .C-13: 
F6:=FCR,C-13:  F7:=FCR-1,C-13: 
CASE  COIR  OF 
0:  BEGIN 

IF  C=3=0)  AND  (F2=255)  THEN  BEGIN  C:=C*i: 

coir:=i; 

END 

ELSE  IF  (F2=0)  AND  (Fl=255)  THEN  BEGIN  R:*R-1 

store;  c:=c*i:  cdir:=o;  end 
else  if  (f1=0)  and  (f0=255)  then  begin  r:=r-1 

cdir:=o;  end 
else  1=  cf0=0)  and  (f7=255)  then  begin  c:=c-1 

store:  r:=r-i:  coir:=3:  end 
else  if  cf7=0)  and  (f6=255)  then  begin  c:=c-1 

cdir:=3s  end 
else  if  (f6=0)  ano  cf5=255)  then  begin  r:=r*1 

store:  C:=c-i:  coir:=3:  end 
ELSE  BEGIN  r:=r*i;  C0IR:=2:  end: 
eno; 
l:  BEGIN 

IF  (F5=0)  ANO  (F4=255)  THEN  BEGIN  R:=R»l: 

CDIR:=2:  ENO 

else  if  (f4=0)  and  cf3=255)  then  begin  c:=c*1 

store:  r:=r*i:  cdir:=i:  end 
else  if  (f3=0)  and  (f2=255)  then  begin  c:=c-h 

cdir:=i:  eno 
else  if  (f2=0)  and  (fl=255)  then  begin  r:*r-1 

store;  c:=c*i;  cdir:=o;  end 
else  if  (f1=0)  and  cf0=255)  then  begin  r:=r-1 

cdir:=o;  eno 
else  if  (f0=o)  and  (f7=255)  then  begin  c:=c-1 

store:  r:=r-i;  coir:=o:  eno 
ELSE  BEGIN  C:=C-i:  COIR:=3:  End; 
eno; 

2:  BEGIN 
•  IF  <F7=0)  AND  (F6=255)  THEN  BEGIN  C:=C-l: 

CDIR:=3:  END 
ELSE  IF  (F6=0)  AND  (F5=255)  THEN  BEGIN  R:=R*1 

store:  c:=c-i:  coir:=3;  end 

ELSE  IF  (F5=0)  AND  CF4=255)  THEN  BEGIN  R:=R*1 
C0IR:=2:  ENO 

else  if  (f4=0)  ano  cf3=255)  then  begin  c:=c*1 

store:  r:=r*i:  coir:=i:  end 
else  if  (f3=0)  and  (f2=255)  then  begin  c:=c+1 

COIR:=i:  END 
ELSE  IF  (F2=0)  AND  (Fl=255)  THEN  BEGIN  R:*R-1 

store:  c:=c*i:  cdir:=i;  end 
else  begin  r:*r-i;  cdir:=o:  eno: 
end: 
3:  BEGIN 

IF  CF1=0)  AND  (FO=255)  THEN  BEGIN  R:=R-i: 

COIR:=0:  END 
ELSE  IF  CF0=0)  AND  (F7=255)  THEN  BEGIN  C:=C-1 
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Source  Listing 


10-Dec-198<.     16:37:51 
10-Dec-198<V     16:37:39 


VAX-1 

_DRA0 


store;  r:=r-i:  coir:=o:  eno 

ELSE  Ic  CF7=0)  AND  (F6=255)  THEN  BEGIN  C:=C-l; 

CDIR:=3:  ENO 
ELSE  IF  (F6=0)  AND  (F5=255)  THEN  BEGIN  R:=«+l; 

store:  c:=C-i:  cdir:=2:  end 
else  if  cf5=0)  and  (f4=255)  then  begin  r:=r+i: 

coir:=2:  eno 
else  if  cf*,  =  0)  ano  (f3=255)  then  begin  c:=c+1; 

store;  r:=r*i;  cdir:=2;  eno 
else  begin  c:  =  c-h;  coir:=i:  end; 
eno; 
end: 
end; 


procedure  initial; 

BEGIN 

FOR  I  :  =  0  TO  255  DO  BEGIN 
XCI3  :=o; 
rciJ  :=0: 
end: 
i:  =  o: 
end; 


PROCEDURE  FIRST: 
VAR 

n,m  :integer; 
BEGIN 

n  :=  o :  m:  =  0  : 
FOR  C  :=0  TO  200  00  BEGIN 
FOR  R  :=0  TO  63  00  BEGIN 

IF  (FC«,C3=255)  ANO  CN=0)  THEN  BEGIN 

Cl  :=  C:  Rl  :=R;  n  :=  N+l ; 
end; 

IF  (FCR,255-C3=255)  ANO  (M=0)  THEN  BEGIN 
C2  :=  255-C:  R2  :=  R;  M  :=M*i; 

end; 
end; 
eno; 
end: 


PROCE 

BEGIN 

NEG 

IF 

F 

E 

ENO 

ELS 

IF 

F 

dure  rotate: 
:  =  o: 

R1  =  R2  THEN   BEGIN  THETA:=0.0: 
OR  l:=0  TO  M  00  BEGIN 

xci3:=c0lci3;  tci3:=ro«ci3: 
no: 

e  theta:=abscarctan((r2-r1)/(c2-c1))); 
cthetaoo)  ano  cr2>r1)  then  begin 

OR  l:=0  TO  M  00  BEGIN 
XCIJ:=C0LEI3*C0SCTHETA)-R0WCI3*SINCTHETA) 
TCI3:=C0LCI3*SINCTHETA)*R0UCI3*C0S(THETA) 
IF  TCI3>=63.0  THEN  NEG:=NEG*l: 

no; 
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Source  Listing 


10-Gec-198<. 
10-0ec-198<. 


16:37: 51 
16:37:39 


VA  X-l 

ctrao 


END 

ELSE  IF  CTHETAOO)  ANO  (R2<R1)  THEN  3EGIN 
FOR  l:=0  TO  M  DO  BEGIN 

xcij:=cqlci3sc0sctheta)*r0wci3ssinctheta); 
ycid:=-calci3*sin(theta)*rouci3*c0sctheta); 
if  yci3>=63.0  then  neg:  =  neg«-1; 
end; 

end: 

(*  FILTER  *> 

I:  =  0; 

for  k:=0  to  m  00  begin 

tempx:=xck.3:  tempy:=yc*3: 

IF  I>0  THEN  BEGIN 

FOR  j:=0  TO  H  DO  BEGIN 

if  (xcj3=tekpx)  and  (ycjd=tempy)  then  i:=j: 
end: 
xC!3:=tempx;  ycI3:  =  tempy  :  m:  =  i:  i:=i*i: 

ENO 

ELSE  BEGIN  XCI3:=TEMPX:  yci3:=tempy:  h:=o;  i:=l 
end; 
end: 
end: 


proc 
m: 
S: 
va 
IP 

PROC 
PROC 
YR 
PROC 
PROC 
PROC 
PROC 
PROC 


EDURE  PARAHC  x:OAi;  Y:DAl;  VAR  Z:OAi:  w:OAi: 
integer:  var  zb:real:  var  ze:real;  k:integer; 
real:  var  n:integer;  var  t:da2:  var  cx:oa2: 
R  Cy:da2:  var  fp:real:  iopt:integer : 
ar:integer:  var  ier:integer);  Fortran: 
edure  initt(speed:integer) :  fortran; 
edure  vwindo(xmin:real:  xrange:real:  ymin:real; 
ange:read:  Fortran: 
edure  move ac x: real:  y:read:  Fortran: 
draha(x:real ;  y:real):  Fortran: 
anchocicharrlnteger) :  fortran; 
finittcii:integer;  I2:integer) 


EDURE 
EDURE 
EDURE 
EDURE 


FORTRAN; 
DASHACXIREAL :  Y:REAL:  L:lNTEGER):  FORTRAN; 


function   DERivc*R=F   t:da2:    n:integer:    Cx:oa2; 

nm:integer;    nu:integer;   arg:real: 

l:integer)    :real:   Fortran: 
proceoure   splcoef; 
begin 

l:=5;    lump:=o: 

WHILE  CK  =  N-5)  AND  (LUMP=0)  DO  BEGIN 

IF(CYCI-13>CYCI3)  AND  ( C YC I ♦ 1 3>C YC I  3 )  AND 
(CYCI*23>CYCI-*13)  THEN 
lump:=i : 
l:»I*l : 
eno: 
if  lu*p=1  then  begin 

flag_begin:=o:  ak:=i;  flag_End:=o;  i:=5: 
uhilI  K  =  N-5  do  begin 

IF  FLAG_BEGIN=0  THEN  begin 

if    ccyci-13>cyci3)   and   (  c  y  c  i  ♦  1  3>c  y  ci  3  )   and 

(cyci*23>cyci*13)   then  begin 

3egcak3:=tci3:   fl ag_segi n : = 1 :   segncak3:=i: 
eno; 
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Source  Listing 


10-Dec-198<.  16:37:51 
10-uec-198<.  16:37:39 


VAX-  1 
CTR  C  : 


END 

else  if  flag_begin=1  then  begin 

if  <cyci-13>=cyci3)  and  ( c y c i 3>  =  c y c i ♦ 1 3 )  and 
<cyci*23>=cyci+13>  then  begin 
encak3:  =  tci-h3:  flag_end:*i;  ennc ak 3 : = i ♦ 1 : 
end: 
end; 
if  cflag_begin  =  i  )  ano  cflag_end=d  then  begin 

flag_bIgin:=o:  flag_end:=o;  ak:=ak*i;  1 :  =  1  - 1 : 
end: 
i:=i+i: 
end: 

END 

ELSE  IF  LUHP=0  THEN  BEGIN 

flag_begin:=o:  flag_end:=o;  i:  =  s:  ak:=i: 
WHIlI  K=N-5  do  BEGIN 

if  flag_begin=0  then  begin 

if(cyci-13>=cyci3>  and  (  c t c i ♦  1  3 >c y c i  3 )  and 

ccyci+23>cyci3)  then  begin 

3egcak3:=tci3:  flag.begin : = 1 :  8egnc  ax.  3  :  =  i  ; 
end: 

END 

ELSE  IF  FLAG_BEGIN=1  THEN  BEGIN 

IF(CYCI-13>  =  CYCI3)  AND  (  < C Y C  I  3< =C YC I ♦ 1 3 ♦ T OL E ) 
ANO  (CYCI3>=CYCI»13-TOL£))  AND 
(CCYCI3<  =  CYCI-»23*TOLE)  AND 
(CYCI3>=CYCI*23*TOLE))  THEN' BEGIN 
ENCAK3:=TCI3;  FLAG_ENO:=i:  ENNCAK3:=i: 

end: 
end: 
if  <flac_begin=1>  ano  cflag_end=1)  then  3egin 

flag_&EGIN:=o :  fla&_eno:=o;  ak:=ak*i:  I:=I-i: 
end; 
i:=I*i: 
end; 

if  flag_begin=i  then  begin 
i:=begncak3: 

while  ck  =  n-5)  and  cflag_end=0)  do  begin 
if  ccyci-13>=cyci3)  ano 
(ccyci  j<  =  cyci-»13*tole) 
and  (cyci3>=cyci*13-tole))  then  begin 
encak.3:  =  tci3:  enncak3:  =  i  :  flag_eno:  =  l: 
end; 
I  :  =  I  + 1 ; 
end: 
end: 
end: 
end: 

procedure  ALUMP; 
BEGIN 
AK  :  =  l ; 
WHILE  (ENNCAK3)>0  00  BEGIN 

AREAC An3:=0 ;  CYMAXCAK3:  =  0.0  ;  C YMI NC AK 3 :  =  1  000 . 0 : 
FOR  I  :  =  (8EGN£AR3)  TO  (ENNCAK3)  DO  BEGIN 
IF  CYMAXCAK3  <=  CYCI3  THEN  BEGIN 
CYMAXCAK3:=CYCI3:  MAXCYCAK3:=TCI3 


3-35 


10-Dec-196<.     16:37:51  VAX-1 

Source    Listing  10-Dec-138-     16:37:39  _D.RAC 

ENO  ; 

if  cymincak3  >=  cyci3  then 
c y m i n c a k3 : = c t c i  3 : 
end; 

cyl:=cycbegncak33-cymincak3: 
for  i  :  =  (begncak3)  to  (enncalo-1)  do  begin 

dct:=ctci*i:-cycid: 

ocx:=cxci-h3-cxci3: 

areacak3:=areacak3*cy1*ocx*ocy*ocx/2.0: 

writelnc -area=*, areacak3.  *  i  =  *,i,*  ak=*,ak): 

cyi:=cyi*dcy; 

ENO; 

ak:  =  ak-h: 

end: 
end: 
procedure  pqsinx; 

BEGIN 

l:=o;  met:=o: 

WHILE  (K  =  M-1)  AND  (MET=0)  DO  BEGIN 
IF  TEMPX=ZCI3  THEN  BEGIN 

tempy:=xcI3:  met:=i: 
end; 
I :=I*1 : 

END; 

end; 

(*******  3**  *£***- *****»*£*<:* 3  *******) 
(*  MAIN  PROGRAM  *) 
BEGIN 

WRITELNC 'INPUT  CUT  OR  CONLINE  FILENAME'); 
READLN(NAME) ; 

OPEN  (INFILE. NAME, HISTORY  :=  OLDi 
ACCESS_METHOD  :=SEQUENTIAL, 

RECORC.LENGTH  : =256 , RECORD _TY PE  :=FIXE0); 
RESETCINFILE): 
R  :  =  0: 

WHILE  NOT  EOF  (INFILE)  DO 
BEGIN 

«EAO  CINFILE.IMAGECR3) I 
FOR  C  :=  0  TO  255  00 

FCR+1,C*13  :=  IMAGECR.C3: 
R  :=  R  +  l ; 

end: 

CLOSc(INFlLE) ; 

first: 

R  :=Ri:  c  :*Ci:  CDIR  :=i: 

initial: 

while  ccoc2)  do  begin 

store: 

cmove: 
end: 
rotate: 

FOR  l:=0  TO  M  00  BEGIN 

if   negoo   then     y  ci  3  :  =  yc  i  3-2  0.  0  ; 

eno: 

m:=m*i;    k:=2:    S:=0.i;    iopt:=o:    ipar:=o: 

for    l:=0    TO    M    00 
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Source  Listing 


10-Dec-l 9f 
10-Dec-19l 


16:37:51 
16:37:39 


VA  X-l 
_QSfiC 


rfCi3:=i.o: 
ans:  =  i  : 
while  ans=1  00  begin 

uritelncold  s=*,s,*new  s='); 

READ(S); 

uritelnc'old  k  =  *  iK« *k=*)  : 

READ(K): 
PARAMCx,Y,Z,W,M,ZB,ZE,K.,S,N,T,CX,CT,FP,IOPT, 

IPAR.IER); 
WRITELNC*  S=*,S.*  IER=*,IER,*  M=',M,'  H=*,N, 

•  CXCN3=',CXCNJ,  *  C»CN-43=*,CTCN-^3): 
WRITELNC  PRINTING   YE5=1  PLOT  X  - Y  =  2  ,C X  ,  C Y  =  3 

'X-Y-CX-C Y=4 ", *  N0=5"); 
READ(ANSl) : 
IF  ANS1=1  THEN  BEGIN 

FOR  l:=l  TO  N  DO 

WRITELN(*CX=*,CXCI3,*    CY»",CYEI3, 
T=',TCI3); 

ENO 
ELSE  IF  ANS1=2  THEN  BEGIN 

INITT(960): 

VWlNOO(0.0,ZCr-- 13, 0.0, 256.0): 

MOVEACZC03.XC03) : 

FOR  l:=0  TO  M-l  00 
0RAWA(ZEI3,XCI3) : 

M0VEA(ZC03,tr03): 

FOR  l:=0  TO  M-l  DO 

dashaczci3,t:i3,2): 

FINITTCOtO) ; 

ENO 

ELSE  IF  ANS1=4  THEN  BEGIN 

INITTC960) ; 

VWlNDOC0.O,ZrM-13,0.0,256.0): 

MovEA(zco3,xno3) : 

FOR  l:*0  TO  M-l  DO 

DRAUA(ZCI3,XrI3) ; 
MOVEACZC03.YC03) : 
FOR  l:=0  TO  M-l  DO 

0RAWA(ZCID,YCI3): 
M0VEA(TC43,CXC43): 
FOR  l:«*  TO  N-4  DO 

DASHACTCI3,CXCI3,2): 
FOR  I:=4  TO  N-4  DO  BEGIN 

M0VEA(TCI3-1.5,CXCI3-1.S); 

anchoc  in ) : 
Enc: 

H0VEACTr_43,CYC43): 
FOR  l:=*  TO  N-4  00 

0ASHA(TCI3,CTCI3,2): 
FOR  l:=*  TO  N-4  00  BEGIN 

M0VEA(TCID-1.5,CYCI3-1.5): 

ANCHO(lll): 

eno: 

Dn,In:  =  o.o:  dmax:-=256.  o  :  tempx:  =  dmin; 

WHILE  TEMPX<=ZCM-1 J  DC  BEGIN 
MOVE A( TEMPX.DMIN); 
DR AWA( TEMPX.OMIN ) : 
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Source    Listing 


10-D<?c-198<.    16:37:51 
10-Dec-19B<.     16:37:39 


VAX-1 

_DT?AC 


ORAUACTEMPx ,DMAX ) : 

TEMPX :=tempx*io.o ; 
end: 

dmin:=o.o:  Dmax:=zcm-13:  tempx:=dmin; 
while  tempx<=256.0  do  begin 

hoveacdmin, tempx)  ; 

drawacdmin, tempx)  ; 

drawacohax,  tempx)  : 

tempx :=tempx*io.o ; 
end: 

FINITTCO.O): 

END 

ELSE  IF  ANS1=3  THEN  BEGIN 
INITTC960): 

VWI NO  0(0. 0,TCN3t0. 0.256.0): 
H0VEA(TC3,CXC43)  : 
FOR  l:*>4  TO  H-U     00 

DRAWA(TCI3,CXCI3); 
FOR  l:=4  TO  N-4  DO  BEGIN 

M0VEA(TCI3-1.5,CXCID-1.5); 

ANCHO(iii) : 

end: 

MOVEACTCAJ.CYC-iD): 

FOR  i:=4  TO  N-*.  00 

DASHA(TCI3,CYCI3,2): 

FOR  l:=A  TO  N-4  00  BEGIN 

H0VEA(TCI3-1.5,CYCID-1.5): 

ANCHO(iii): 

end; 
finittco.o): 
end: 
(*  important  fortran  declear  from  1  * ) 

nki:=n-k-i:  l:=k*i:  nEnd:=n-k;  ji:=o: 

FOR  i:=L  TO  NENO  DO  BEGIN 
ARG:=TCI3: 

hhil=(arg>=tcl*13)   and  (l<nm)  do 

l:  =  l*i  : 
spx1cj13:=  0erivct,n,cx,nk1,0,arg,l): 

SPT1CJ13:=  DERIV(T,N,CY,NKl,OtARG,L): 

Jl:  =  Jl-»i; 

end; 

j:  =0  :  l:=k  +  i ; 

for  i:=l  to  neno  co  begin 

ARG:=TCI i; 

while  (arg>=ttl*13)  and  (l<nk1)  do 

l:=l*i: 
tempx:=tci*13-arg: 

if  tehpx>=7.0  then  begin 
while  arg<tci*13  do  begin 

SPXCJ3  :=  DERI V(T,N,CX,NK1,0,ARG,L): 
SPYCJ3  :=  DERIV(T,N,CY,NK.1,0,ARG,L): 

arg:=arg*2.0;  j:=j*1: 

end: 

END 

ELSE  BEGIN 

SPXCJ3  :=  DERIV(T,U,CX,NK1,0,ARG,L); 

SPYCJ3  :=  DERIV(T,N,CY,NK1,0,ARG,L) : 
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Sourc  e  Listing 


10-0ec- 
10-0ec 


19B<. 
193^ 


16:37: 51 
16:37: 39 


VAX-] 

_DR  AC 


PLOTING 


Y£S  =  1 


,5) 


J  :  =  J  ♦  1  : 

end: 
end: 
writelnc'do  tou  want 

*print=2  compare=3 
readcans2) : 
if  ans2=1  then  begin 

INITTC960): 

VUINDOCO.0.2  56. 0,0. 0.256.0): 

H0VEA(XC03,TC03) : 

FOR  l:=0  TO  M-l  00 

DRAWA( XCIl.YClJ}; 
FOR  l:  =  0  TO  Jl-1  DO  BEGIN 

HOVEACSPX1CI3-1.5.5PY1CI3-1 

ANCHGClll): 

end: 
m0veacspxc03.sptc0]): 

FOR  l:=l  TO  J-l  00 

0ASHA(SPXCI3.SPYCI3,2): 
FINITTCO ,0): 
END 

ELSE  IP  ANS2=2  THEN  BEGIN 
FOR  l:=0  TO  Jl-1  DO 

WRITELM(*SPX1=',SPX1CI3, 
SPY1=-,5PY1CI3): 
END 

ELSE  IF  ANS2=3  THEN  BEGIN 
WRITELNC  *TOL=*)  ; 
READCTOLE): 
SPLCOEF; 

aluhp; 

ak:=i:  ra:=xcm-i J-XC03: 

WHILE  EN:aK.3>0  DO  BEGIN 
TEMPx:-BEGCAK.D  :  POSINX: 
STARCAK3:=(CTEMPY-XC03)-RA/2)/ 

tempx:=enc AK3:  posinx: 

TERCAK3:=<CTSMPY-XC0D)-RA/2)/R 

tempx : =maxc YCAX3  :  posinx; 

PEAK.CAK;j:  =  (CTEHPr-XC0j)-RA/2)/ 
AREACAK.J  :  =  (AREAC  AHD)/(R  A**2  )  : 
WRITELNC -EEGIN  =  ',STARCAK3f  •  EN 
TERCAK.J,'  TOTAL=*,RA, 
'  AREA=*  .AREAC AK3, 
"  PE AK  =  '  .PEAKCAK])  ; 
AK :=AK*1 ; 

eno: 
eno; 
writelnc'do  tou  want  run  again  yes 

*  N0=2*); 

reaccanS); 
if   ans=1   then   eegin 
writelnciqpt  =  -)  : 
readciopt): 
end: 
end:    (*    WHILE    *o 
END.    ' 


ra  : 

a  : 
r  a: 

0=', 


=  1-, 
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